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I.  INTRODUCTION 

Various  methods  of  four-dimensional  data  assimilation  have  been  in 
practice  both  operationally  and  in  a  data  archival  mode  over  the  past  ten 
years.  Most  current  methodologies  use  a  global  objective  analysis 
procedure,  an  initialization  program  to  remove  gravity  modes  and  their 
tendencies,  and  a  short-term  (usually  6h)  global  forecast  model  to  step 
ahead  to  the  next  time  period.  An  example  of  an  operational  global  data 
assimilation  system  is  that  used  at  the  National  Meteorological  Center  (NMC) 
(Dey  and  Morone,  1985),  while  the  European  Centre  for  Medium  Range  Weather 
Forecasts  (ECHWF)  has  used  its  data  assimilation  procedure  for  the 
formulation  of  the  FGGE  III-B  analysis  archives  (Bengtsson  et  al. ,  1982)  as 
well  as  for  operational  purposes  (Lorenc,  1981). 

Until  recently,  the  focus  of  attention  in  the  development  and  use  of 
four-dimensional  data  assimilation  systems  has  been  on  mass  and  motion 
fields.  Although  the  NMC  objective  analysis  procedure,  the  optimum 
interpolation  method,  does  use  analyzed  winds  and  temperatures  in  its 
special  univariate  horizontal  and  vertical  correlation  functions  for 
relative  humidity,  the  observational  data  used  in  the  analysis  are  limited 
to  rawinsonde  and  surface  weather  observation  measured  humidities. 
Operationally,  ECHWF  has  analyzed  integrated  water  vapor  mixing  ratio  of 
layers  using  measured  temperature  and  dewpoint  from  rawinsonde  reports  and 
relative  humidity  inferred  from  four  layers  below  300  mb  from  surface 
weather  observations  (Tibaldi,  1982).  Recently,  the  use  of  layer 
precipitable  water  from  satellite  soundings  has  been  studied  at  ECHWF 
(Lonnborg,  private  communication)  but  as  yet  it  has  not  been  included  in  its 
operational  optimum  interpolation  procedure  for  moisture.  Therefore,  the 
large  regions  of  the  globe  where  conventional  observations  of  humidity  are 
sparse  or  non-existent  are  not  updated  in  the  analysis  phase  of  the 
assimilation,  and  in  these  regions  the  forecast  model  is  allowed  to  create 
its  own  moisture  climatology. 

An  accurate  depiction  of  cloud  distribution  is  of  great  importance  in 
Air  Force  operations,  in  both  the  analysis  and  forecast  modes.  Thus,  it  is 
essential  that  some  kind  of  updating  of  the  forecast  humidity  field  take 
place  in  data  void  regions  to  ensure  the  best  possible  representation  of 
clouds  in  those  regions  from  both  the  present  analysis  and  forecast 
humidities.  Obviously,  the  expense  of  attempting  to  cover  data 


void  regions  adequately  with  conventional  observations  on  a  routine  basis 
would  be  prohibitive,  so  other  forms  of  humidity  information  must  be  sought 
for  such  regions.  Previously  mentioned  satellite  profiles  of  layer 
precipitable  water  offer  some  hope  if  they  are  recovered  from  satellite 
radiance  measurements  with  sufficient  accuracy.  Another  form  of  global 
moisture  information  that  may  prove  useful  in  updating  forecast  humidities 
is  relative  humidities  inferred  from  cloud  cover  values  of  automated 
nephanalyses .  One  source  of  cloud  analyses  is  the  3DNEPH  cloud  analysis 
procedure  of  the  Air  Force  Global  Weather  Central  (AFGWC)  (Fye,  1978). 

This  report  focuses  on  the  use  of  such  alternative  forms  of  humidity 
information  to  supplement  conventional  humidity  observations  in  preparing 
global  analyses  of  relative  humidity  in  a  data  assimilation  environment. 
Humidities  derived  from  satellite  soundings  and  inferred  from  surface 
weather  observations  (SWO)  and  AFGWC  cloud  analyses  (3DN)  are  studied  for 
possible  use  in  the  AFGL  Statistical  Analysis  Program  (ASAP),  and  data 
assimilation  sequences  are  run  with  and  without  the  more  promising  forms  of 
alternative  humidity  information.  Conclusions  are  then  drawn  from  the 
results  as  to  the  effect  on  both  humidity  analysis  climatologies  and 
humidity  forecasts  of  using  such  supplemental  humidity  information. 

II.  DESCRIPTION  OF  THE  GLOBAL  AFGL  DATA  ASSIMILATION  SYSTEM  (AGDAS) 


The  triad  of  objective  analysis,  initialization,  and  global  forecast 
model  programs  used  in  the  global  AFGL  data  assimilation  system  consists  of 
the  ASAP  optimum  interpolation  (01)  analysis,  the  Machenhauer  non-linear 
normal  mode  initialization  (NMI)  (Ballish,  1980)  using  four  vertical  modes 
and  two  iterations,  and  the  AFGL  global  spectral  model  (GSM)  (Brenner  et 
al. .  1982).  The  AFGL  GSM  used  in  this  study  was  modified  in  the  structure 
of  its  numerics  from  the  NMC  spectral  model  (Sela,  1980),  but  changed  in  the 
physics  only  to  allow  for  a  different  distribution  of  vertical  layers.  Two 
later  upgrades  of  the  physics  parameterizations  (NMC  1983,  NMC  1985)  have 
since  been  carried  out  in  the  NMC  operational  GSM.  A  30  wave  rhomboidal 
truncation  on  12  a  layers  is  used  in  the  NMI  and  GSM,  and  the  analysis  is 
performed  on  a  layers  in  the  ASAP  to  avoid  vertical  interpolation  of 
fields  between  pressure  levels  and  a  layers.  The  remainder  of  the 
description  in  this  section  concerns  the  01  procedure  ASAP. 

Originally,  ASAP  was  developed  directly  from  the  1979  multivariace  01 
procedure  as  described  by  Bergman  (1979);  in  fact,  the  codes  were  developed 


using  only  Bergman's  article  and  the  companion  paper  by  McPherson  et  al.. 
(1979)  as  references.  With  some  assistance  from  MMC  personnel  (Morone, 

1983,  private  communication),  a  multivariate  temperature  (T)  and  zonal  (u) 
and  meridional  (v)  wind  component  analysis  code  on  a  layers  and  a  separate 
surface  pressure  01  analysis  code  were  formulated.  The  surface  pressure 
analysis  was  later  dropped  and  the  T,  u,  v,  a  layer  analysis  was  replaced 
by  a  height  (Z),  u,  v,  analysis  on  a  layers.  The  equations  for  computing 
the  corrections  for  an  analysis  grid  point  based  on  weighted  sums  of 
surrounding  observation-minus-first  guess  residuals  for  the  normalized 
analysis  error,  and  for  the  computation  of  the  horizontal  and  vertical 
correlation  functions  from  Bergman  are  used  without  modification.  Still 
later,  cycling  of  the  analysis  error  and  the  great  circle  distance  method 
for  correlation  functions  equatorward  of  70*  latitude  were  included  as 
described  by  Dey  and  Morone  (1985)  without  changing  the  Bergman  formulation 
(including  map  factor)  for  latitudes  |+|  >  70*.  The  following  paragraphs 
give  a  detailed  description  of  the  four  codes  that  make  up  ASAP. 

A.  The  Pre-Analysis  Program  (PRAMAL) 

PRAMAL  converts  spectral  coefficients  of  absolute  vorticity  (n) , 
divergence  (D),  T,  specific  humidity  (q,  carried  only  in  the  lowest  seven 
0  layers),  and  surface  pressure  pA  from  the  previous  forecast  (first 
guess)  to  grid  point  values  on  both  the  analysis  grid  (61  equally  spaced 
longitudes  by  62  Gaussian  latitudes)  and  the  interpolation  grid  (1*  latitude 
x  1*  longitude).  Thus,  values  of  T,  u,  v,  relative  humidity  (r),  and  pA 
are  located  on  the  respective  grids  on  12  a  layers  for  T,  u,  v,  and  seven 
a  layers  for  r,  representing  the  first  guess.  In  PRAMAL,  each  grid  value 
of  T,  q,  and  pA  is  used  to  compute  the  corresponding  grid  value  of  r  using 
the  relations 
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where  the  subscript  &  represents  the  index  of  the  a  layer.  In  addition, 

PRAKAL  also  converts  the  spectral  coefficients  of  height  analysis  error 

(a_)(12  0  layers)  and  relative  humidity  analysis  error  a  (seven 
z  r 

o  layers)  from  the  previous  analysis  to  61x62  grid  values.  In  cases  where 
the  present  cycle  is  the  first  in  a  series,  the  analysis  error  conversion  is 
not  performed  in  PRAHAL,  since  no  previous  analysis  in  that  series  exists. 

B.  The  Mass-Motion  Analysis  Program  (ASAPOI) 

The  upper  air  analysis  of  Z,  u,  v  is  multivariate  and  is  carried  out  at 
the  sigma  layer  surfaces  of  the  GSM  forecast.  Fig.  1  shows  the  relative 
positions  of  the  sigma  layers  (0^)  and  sigma  levels  (o^)  currently  used 
in  the  AGDAS .  Horizontal  (w)  and  vertical  (u)  correlation  functions  for 
the  autocorrelations  and  cross-correlations  involving  Z,  u,  v  used  are  those 
given  by  Bergman  (1979)  for  latitudes  $  such  that  |$|>70*,  and  by  Dey 
and  Morone  (1985)  for  |t|<70*.  The  respective  distance  formulations  from 
the  two  articles,  as  trail  as  the  values  for  constants  from  Dey  and  Morone, 
are  used.  ASAPOI  also  uses  the  rawinsonde,  pibal,  and  single-level  wind 
observation  error  standard  deviation  values  from  Dey  and  Morone,  but  uses 
standard  deviation  values  for  satellite  height  observation  errors  derived 
from  a  collocation  study  with  rawinsondes  to  be  discussed  later  in  this 
report.  Values  for  k^,  the  constant  in  the  horizontal  autocorrelation 
function  for  Z,  are  the  same  as  those  specified  by  Dey  and  Morone,  but  are 
used  at  0  layers  1-9,  0  layers  10  and  11,  and  0  layer  12  respectively. 

Because  the  analysis  is  performed  on  the  0  layers  rather  than  on 
pressure  surfaces,  representations  of  the  various  types  of  observations  used 
in  the  analysis  must  be  obtained  at  the  pressure  levels  of  the  0  layers  at 
each  observation  site.  Interpolation  procedures  used  for  this  purpose  are 
discussed  in  detail  in  the  following  section.  ASAPOI  inputs  four  types  of 
observations  from  the  FGGE  II-B  data  sets:  (1)  Type  1  observations,  which 
include  rawinsondes,  pilot  balloons  (pibals),  dropsondes  and  other  direct 
profiling  instruments,  (2)  Type  4  observations,  which  are  the  satellite 
soundings  of  layer  mean  temperature  and  layer  precipitable  water,  (3)  Type  2 
observations,  aircraft  observations,  and  (4)  Type  6,  single  level  cloud 
track  wind  observations.  Each  data  type  is  ordered  by  latitude  from  south 
to  north  prior  to  its  use  in  ASAPOI,  in  order  to  facilitate  the  location  of 


the  proper  latitudes  of  1*  x  1*  first  guess  values  from  which  to  interpolate 
bilinearly  to  each  observation  site  as  it  is  processed.  This  horizontal 
interpolation  to  the  observation  location  is  performed  for  the  a  layer 
variables  T,  u,  v,  r,  and  for  p#.  Then  a  layer  T  is  converted  to  a 
level  Z  using 


E*  +  (R/g)T  ln(o.  ./a.) 
l-l  °t  1-1  1 
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where  Z*  3  model  terrain  height  interpolated  from  the  1°  x  1* 

grid,  anS  a  >0.01  for  the  purposes  of  establishing  a  value  of  Z  for  the 
uppermost  level.  These  values  of  Z,  u,  v,  r  are  subtracted  from  the 
corresponding  estimates  of  observation  values  of  Z,  u,  v,  r  to  form  the 
residuals  used  in  the  analysis.  The  residuals  are  stored  sequentially  in  a 
data  list  in  the  order  they  are  processed,  from  south  to  north  for  each  data 
type  in  turn. 

ASAPOI  then  calls  a  quality  control  procedure  that  performs  a  "gross" 
check  and  "buddy"  check  of  the  Z,  u,  v  residuals.  Residuals  are  sorted  by 
location  into  buddy  check  boxes  that  are  10*  in  latitude  extent  and 
approximately  equal  in  longitude  extent,  with  the  number  of  boxes  for  each 
latitude  band  equal  to  NBOXRS  =  36*  cos  +,  where  t  is  the  central 
latitude  of  each  band.  For  boxes  whose  central  latitude  lies  poleward  of 
70*  latitude,  all  wind  component  residuals  are  converted  from 
latitude-longitude  (♦,  X)  orientation  to  Cartesian  (x,  y)  orientation 
using 


d“x,y  *  -%.X  *in  x  -  Wx  c“ 
d,x.y  -  Wx  C°*  X  ’  “Vx  sin  X 


X 


(6) 


in  which  du,  dv  represent  the  wind  component  residuals,  X  the  observation 
longitude,  and  1=1  for  Northern  Hemisphere  and  -1  for  Southern 
Hemisphere.  Next,  the  gross  check  is  performed  for  each  box  by  first 
computing  the  forecast  error  standard  deviation  for  each  individual  residual 
at  its  respective  level  (Z)  or  layer  (u,  v)  by  vertically  interpolating  the 
values  in  Table  1  to  the  a  pressure.  If  a  residual  exceeds  three  times 
the  size  of  the  forecast  error  standard  deviation  value  interpolated  from 
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Table  1  it  is  removed.  If  either  wind  component  fails  this  gross  check, 
both  are  removed. 

The  buddy  check  is  then  performed  univariately  within  each  box  at 
each  level  by  computing  total  autocorrelations  for  each  respective  residual, 
and  determining  if  the  criteria 

|d2i  -  dZjl  <  (a  -  bp  “)4  J 

|dUi  -  dUj  |  <  (a  -  bp  ““)4  *  (7) 

|dvA  -  dVj|  <  (a  -  bp  £^)4  J 

are  met.  In  these  expressions,  i  and  j  are  the  indices  of  two  observations 
being  compared,  4^ ,  4^  ,  4^  ,  are  the  forecast  error  standard  deviation  values 

computed  from  Table  1  for  the  respective  variables  at  the  location  i, 

Pjj  are  the  spatial  autocorrelations  and  as  3.0,  b  *  1.5.  If  the 
criteria  are  met  for  the  comparison  of  residuals  and  i  and  j  and  p^  > 

0.75,  the  observation  considered  to  be  of  higher  quality  is  assigned  a  keep 
flag;  if  they  are  of  equal  quality,  both  are  assigned  keep  flags.  For  wind 
residuals,  p^  for  u  and  v  must  be  >  0.75.  If  the  criterion  is  not  met, 
the  lesser  quality  observation  is  assigned  a  toss  flag,  and  both  get  a  toss 
flag  if  they  are  of  equal  quality.  Separate  keep  and  toss  flags  for  wind 
components  are  not  maintained — a  wind  keep  flag  is  assigned  if  p^  for 
both  u  and  v  is  >  0.75,  and  a  wind  toss  flag  is  assigned  if  either  u  or  v 
fails  the  buddy  check  criterion.  The  0.75  factor  is  used  to  insure  that 
only  more  spatially  correlated  observations  corroborate  each  other. 

Quality  level  in  the  buddy  check  is  defined  by  observation  type  and 
whether  or  not  a  vertical  consistency  check  was  conducted  for  that 
observation.  Only  FGGE  II-B  observations  with  quality  control  indices  of  0 
(no  check  made)  and  1  (observation  found  correct  in  vertical  consistency 
check)  are  used  in  the  analysis,  with  the  latter  given  the  higher  quality 
level.  The  quality  control  indicator  (QCI)  assigned  to  a  residual  for  the 
purposes  of  determining  comparative  quality  is  set  as  follows:  QC  >  0  for 
vertically  checked  data  and  =  1  for  vertically  unchecked  observations.  Then 


for  Typo  1  data,  QCI  *  QC+1,  while  QCI  for  Types  2,  4,  and  6  are  given  by 
QCI  -  3(QC+3)+ITBI,  where  ITBI  «  1  for  Type  2,  2  for  Type  4,  and  3  for 
Type  6.  For  example,  a  QC  value  of  0  will  correspond  to  QCI  -  1  for  Type  1, 
10  for  Type  2,  11  for  Type  4,  and  12  for  Type  6.  Thus,  the  higher  the  value 
of  QCI,  the  poorer  the  quality  assumed  for  the  residual.  In  this  way,  the 
conventional  upper  air  data  (Type  1)  will  never  be  flagged  by  other 
observation  types. 

Once  all  of  the  observation  intercomparisons  are  made  in  a  particular 
box  for  a  particular  level,  each  observation  is  checked  to  see  if  it  has  two 
or  more  keep  flags.  If  so,  all  toss  flags  are  removed  from  the  observation 
since  it  is  to  be  retained  regardless  of  how  many  toss  flags  it  has 
received.  Observation  iterative  rejection  begins  with  summation  of  toss 
flags  assigned  to  each  residual,  and  the  residual  of  each  type  (heights  and 
winds  separately)  with  the  largest  number  of  toss  flags  is  identified.  This 
residual,  along  with  all  others  that  have  within  four  of  that  maximum 
number,  is  eliminated,  and  the  toss  flags  that  they  imposed  on  the  remaining 
residuals  are  removed .  The  remaining  residual  of  each  variable  with  the 
largest  number  of  toss  flags  is  identified;  it  and  others  within  four  toss 
flags  of  it  are  eliminated,  the  imposed  toss  flags  on  other  residuals  are 
removed,  and  so  on.  The  four  flag  check  is  removed  once  the  highest  number 
of  toss  flags  falls  below  four.  This  continues  until  all  residuals  of 
height  and  winds  have  at  most  one  toss  flag,  or  ten  iterations  have  been 
carried  out,  whichever  comes  first.  Remaining  residuals  are  candidates  for 
the  01  analysis.  The  wind  residuals  in  boxes  whose  central  latitude  is 
poleward  of  70*  are  converted  back  to  +,  X  orientation. 

A  number  of  points  (PTS)  is  assigned  to  each  observation,  whether 
sounding  or  single  level.  For  each  o  level  of  a  particular  observation, 
if  a  height  observation  is  present,  PTS  is  added  to  the  sum  for  the 
observation.  For  wind  observations  at  each  a  layer,  one  point  is  added  to 
the  sum,  then  if  the  wind  observation  is  accompanied  by  a  height  observation 
at  the  sigma  level  just  above  it,  an  additional  value  of  PTS  is  added. 

Here,  PTS  has  a  value  of  1.0  for  all  observation  types  except  Type  4,  where 
PTS  -  0.42.  This  lower  factor  represents  the  overall  reliability  (ratio  of 
observation  error  standard  deviation  for  rawinsonde  heights  to  that  computed 
for  satellite  sounding  heights)  of  satellite  soundings  compared  with 
conventional  height  observations.  The  total  number  in  points  for  each 


observation  over  all  levels  is  stored  and  used  later  in  the  selection  of 
observations  lying  around  an  analysis  point.  It  can  be  seen  that  the  more 
information  contained  in  an  observation,  the  larger  its  number  of  points — it 
has  more  information  to  contribute  to  a  particular  location  on  the  analysis 
grid. 

At  this  point,  the  set  of  observation-minus-first  guess  residuals 

contains  checked  residuals  of  Z  on  a  levels  and  u,  v  on  a  layers  where 

they  exist  in  each  observation,  as  well  as  unchecked  residuals  of  r  on  a 

layers.  Observations  have  been  assigned  points  and  are  ordered  within  each 

type  by  latitude.  The  analysis  of  the  61x62  grid  points  then  begins  with 

the  southern-most  Gaussian  latitude  and  the  Greenwich  Meridian,  moving 

eastward.  For  each  analysis  point,  the  following  sequence  of  steps  takes 

place  to  determine  the  NSNDS  (currently  8)  observations  that  will  be  used  to 

update  the  first  guess  values  at  all  a  layers .  First ,  the  Type  1 

observations  are  considered,  by  computing  the  horizontal  height 

zz 

autocorrelation  p.  between  each  observation  i  and  the  grid  point  g.  Before 

HSNDS  observations  with  p^  >0.1  are  collected,  any  observations  failing  this 

criterion  are  no  longer  considered.  When  HSNDS  such  observations  are  found, 

ZZ 

each  successive  observation's  value  of  p.^  is  multiplied  by  PTS,  and  this  product 

is  compared  with  the  same  product  for  the  NSNDS  observation,  where  observations 

1 ,  2 ,  ....  HSNDS  already  collected  are  ordered  from  largest  to  smallest  value 

ZZ 

of  this  product.  If  a  subsequent  observation's  value  of  p.  *PTS  is  larger  than 

xg 

that  for  the  last  observation  in  the  collection,  the  latter  is  excluded  from 

ZZ  ZZ 

consideration.  Making  the  test  for  inclusion  on  p,  rather  than  on  p.  *PTS 

ig  ig 

for  the  first  HSNDS  observations  avoids  making  the  amount  of  information  in 
an  observation  a  criterion  for  its  inclusion  for  analysis  points  in  data 
sparse  areas.  However,  in  data  rich  areas  the  combination  of  nearness  to 
the  grid  point  and  amount  of  information  influences  an  observation’s  chances 
for  participation  in  the  analysis.  This  procedure  continues  until  all  Type 
1  observations  have  been  considered.  If  less  than  NSNDS  observations  have 
been  found  from  the  complete  set  of  Type  1  observations,  the  entire  set  of 
Type  4  observations  is  compared  against  the  accepted  Type  l's  in  the  same 
manner.  If  NSNDS  surrounding  observation  sets  still  have  not  been  found, 
the  Type  2  and  Type  6  sets  are  considered  in  order.  Thus  the  search 
continues  through  the  four  data  types  in  this  order,  discontinuing  after 


10 


considering  in  its  entirety  the  first  data  type  set  in  which  the  collection  of 
VSHDS  observations  has  been  completed.  The  entire  set  of  VSHDS  observations  is 
a  candidate  for  participation  in  the  procedure  that  corrects  that  grid  point's 
first  guess  values  at  all  model  layers.  In  the  present  version  of  the 
analysis,  VSHDS  *  8,  but  a  correction  can  take  place  even  if  only  one 

zz 

observation  has  a  >0.1.  If  no  observations  satisfy  this  condition,  no 

correction  is  calculated  for  the  analysis  point  at  any  layer. 

The  purpose  for  assigning  points  to  the  observations  is  to  indicate  how 
much  that  observation  could  contribute  to  correct  the  grid  point  values  at 
all  <t  layers.  For  example,  based  on  the  current  point  system  used,  a 
complete  rawinsonde  observation  for  a  12  layer  model  structure  would  have  36 
points,  5.04  points  for  a  complete  satellite  sounding,  and  1  point  for  an 
aircraft  or  single  level  satellite.  This  means  that  a  complete  rawinsonde 
observation  located  7.14  times  farther  from  the  analysis  point  than  a 
complete  satellite  sounding  or  36  times  farther  than  an  aircraft  or  single 
level  satellite  wind  observation  would  have  equal  likelihood  of  being  chosen 
for  that  grid  point  in  data  rich  regions.  In  data  sparse  regions  where 
there  are  fewer  than  VSHDS  of  all  types  nearby,  the  observation  types  have 
equal  opportunity.  Types  of  data  are  processed  in -the  indicated  order 
because  of  the  perceived  accuracy  of  the  observations  and  the  number  of 
residuals  generally  present  in  each  type.  Conventional  upper  air  sounding 
is  seen  as  the  most  accurate  and  as  having  the  most  residuals.  Satellite 
soundings  are  taken  next  because  of  the  complete  set  of  height  residuals 
they  provide,  then  aircraft  winds  which  are  more  accurate  than  satellite 
winds . 

Once  the  observations  are  collected  for  a  particular  analysis  point, 
the  01  procedure  begins.  For  the  first  analysis  in  the  assimilation 

sequence,  the  forecast  error  standard  deviation  relation  6  *  6  = 

1/2  u  v 
6z[(2kh>  g/fg],  where  k^  is  the  constant  in  the  height  autocorrelation 

function,  g  is  acceleration  due  to  gravity,  and  f  is  the  coriolis 

8 

parameter  at  the  analysis  point,  is  used  to  compute  &u>  &v  poleward 
of  25*  latitude  using  the  4^  values  in  Table  1.  Within  25*  of  the 
equator,  6u,  6^  values  derived  in  this  way  are  blended  with  tropical 
values  following  Dey  and  Morone  (1985).  Then  mandatory  level  values  are 
linearly  interpolated  in  In  p  (natural  logarithm  of  pressure)  to  the  sigma 
layer  pressures.  For  subsequent  analyses  in  an  assimilation  series,  the 
analysis  error  values  for  Z  (<*,)  from  all  layers  of  that  analysis  point 
from  the  previous  analysis  are  read  into  the  code,  having  previously  been 
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converted  from  spectral  form.  Values  for  the  6h  height  forecast  error  growth 
amount  for  each  sigma  layer,  estimated  from  interpolating  the  corresponding 
mandatory  pressure  level  values  from  Dey  and  Korone  (1985)  to  the  sigma 


layers,  are  added  to  to  form  6^.  However,  before  this  addition 


takes  place,  the  az  values  are  checked  to  make  sure  that  they  are  greater 
in  magnitude  than  the  6h  forecast  error  growth  amount,  yet  less  than  the 
climatological  limit  for  height  forecast  error  interpolated  to  the  sigma 
layers  from  the  mandatory  level  Dey  and  Horone  values.  Therefore,  the  6 
values  are  allowed  to  range  between  twice  the  6h  forecast  error  growth  amount 
as  a  lower  limit  and  the  sum  of  the  climatological  limit  and  the  6h  forecast 


growth  amount  as  an  upper  limit.  From  these  values,  the  expression 
1/2 

(2k^)  g/f  was  used  poleward  of  25*  and  blending  with  tropical  values 


was  used  equatorward  of  25*  to  obtain  &u,  6^  values  at  all  a  layers  for  the 


analysis  point.  The  forecast  error  standard  deviation  values  assigned  to  each 


residual  i,  4,  ,  6  ,6  ,  is  the  analysis  a  layer  value  for  the  same 

Zi  Ui  Vi 


variable  at  the  location  of  the  analysis  point  in  the  horizontal  and  at  the 
a  layer  pressure  of  the  analysis  point  closest  to  the  pressure  of  the 
observed  residual  in  the  vertical.  Thus,  local  homogeneity  of  the  forecast 


error  standard  deviation  is  assumed  for  the  horizontal  but  not  for  the 


vertical . 


Next,  the  NLRHS  (presently  10)  individual  a  level  height  or  a  layer 
wind  (both  components  taken  together)  residuals  to  be  used  in  the 
computation  of  a  particular  a  layer's  values  of  corrections  AZ,  Au, 


Av  are  selected  from  the  NSNDS  set  of  observations.  For  each  residual 


Z_  ,  u  ,  v  in  the  complete  set,  all  auto-  and  cross-correlations  with  the 
°i  °i  °i 

2 

sigma  layer  analysis  point  are  computed  and  divided  by  1  +  where  c 


(normalized  observation  error)  is  the  ratio  of  the  observation  error  standard 


deviation  for  that  residual  to  the  forecast  error  standard  deviation  A. 

X  #  #  0  X 


ZZ  7n  Zv 

for  that  residual.  These  divided  correlations  p^  ,  P  ^ .  and  p ^  for  each 


height  residual  Z  and  p?Z,  p?u,  p*fv,  pYZ,  oV1 ,  p^  for  each  wind  residual 

oi  ig  'lg’  ig  ig  rig  ig 


pair  u  ,  v  are  used  as  a  basis  for  the  selection  for  the  NLRHS  height  and 
°i  °i 


wind  residuals  used  to  compute  all  three  of  AZ,  Au,  Av.  It  was  noticed 
that  by  choosing  the  10  closest  residuals  to  the  grid  point  (i.e.,  largest 


zz 

values  of  for  all  residuals)  on  the  basis  of  distance  alone,  many  times 

the  closest  observations  had  the  lowest  actual  correlation  with  the  grid  point. 

2 

By  dividing  the  correlation  by  1  +  ,  the  observation's  reliability  is  used 

as  a  factor  in  its  selection.  Less  reliable  observations,  with  their  larger 

value  of  c^,  will  not  have  as  much  influence  in  the  analysis.  Since  the  NLRHS 

residuals  selected  will  affect  the  computation  of  AZ,  Au,  and  Av,  it  was  thought 

to  be  important  to  consider  the  combined  correlations  of  each  residual  with  grid 

ZZ  *  Zu "  Zv  * 

point  Z,  u,  and  v.  Thus,  the  sum  of  the  absolute  values  of  p,  ,  p,  ,  p. 

*8  *SU2 *  *8 

for  each  height  residual  and  half  the  sum  of  the  absolute  values  of  p ^  , 

nn  iiv  vZ  VU  w 

p.  ,  p,  ,  p.  ,  p.  ,  p.  for  each  pair  of  wind  component  residuals  were 
ig  xg  ig  ig  lg 

compared  against  each  other,  and  the  NLRHS  residuals  with  the  highest  values 

of  these  sums  (all  of  which  had  to  be  >  0.1)  are  those  chosen  for  the 

analysis  at  the  a  layer.  If  fewer  than  NLRHS  residuals  are  found  to  have 

sums  >0.1,  only  those  that  do  are  used.  In  this  tray,  the  number  of 

individual  residuals  (Z,  u,  and  v  separately)  can  be  as  low  as  one  (if  only 

one  Z  residual  is  identified)  and  as  large  as  2  *  NLRHS  (if  HLRHS  wind 

residuals  have  the  largest  sums,  to  the  exclusion  of  any  Z  residuals).  No 

limit  is  imposed  on  the  number  of  height  or  wind  residuals  that  may  be 

selected  from  a  particular  observation.  Because  height  residuals  are 

carried  at  a  levels  and  wind  residuals  are  at  a  layers,  correlations 

with  a  layer  analysis  points  favor  wind  residuals  because  of  their  larger 

ZZ 

value  of  the  vertical  correlation  v.  that  is  used  for  all  residuals.  However, 
this  is  offset  somewhat  by  the  relative  magnitudes  of  c  and  c  at 

a 0  U  .  V 


Values  of  the  ratio  c_/c  shown  in  Table  2  were 

Z  u,v 


most  levels. 

computed  for  three  latitudes  and  12  mandatory  levels  using  the  relationship 


c  /c 
Z  u,v 


.1/2  _ 
(2kh)  8*Z 


(8) 


n 


u,v 


where  values  of  the  observation  error  standard  deviation  for  heights 

(£z)  and  winds  (5  y)  for  rawinsondes  from  Table  5  of  Dey  and 

Morone  (1985)  as  used  in  the  ASAPOI,  are  used.  For  levels  below  200  mb  and 

latitudes  poleward  of  30* ,  the  lower  value  of  cz  reduces  the 

2 

correlations  divided  by  1  +  for  heights  less  than  the  divided  wind 


Table  2.  cz/cu>v  for  Rawinsonde  Observations  for  Three  Latitudes 


30* 

1000 

1.045 

850 

1.1946 

700 

1.1007 

500 

.9566 

400 

.9535 

300 

.9718 

250 

1.1571 

200 

1.2619 

150 

1.4760 

100 

1.6675 

70. 

3.0022 

50 

2.8899 

45* 

60* 

.7391 

.6035 

.8447 

.6897 

.7783 

.6355 

.6764 

.5523 

.6743 

.5505 

.6871 

.5610 

.8182 

.6681 

.8923 

.7286 

1.0437 

.8522 

1.1791 

.9627 

2.1229 

1.7333 

2.0434 

1.6685 

correlations  with  tha  larger  value  of  c  .  The  offsetting  of  these  two 

U  #  V 

factors  led  to  a  good  overall  blend  of  height  and  wind  residuals  used  in  the 
correction  computations  over  the  entire  globe. 

The  coefficient  of  geostrophy  G  used  in  the  horizontal 
height-wind  cross-correlations  equatorward  of  70*  latitude  (see  Oey  and 
Mo rone* s  gqs.  B2-B5)  for  both  observation-analysis  point  correlations  and 
observation-observation  correlations  is  a  blend  of  the  two  proposed  by 
Bergman  (1979)  and  Dey  and  Mo rone  (1985).  It  was  seen  that  the  diminishing 
of  the  coupling  to  zero  between  25*  and  10*  was  important,  yet  the  sudden 
change  of  vs.  latitude  in  the  vicinity  of  25*  was  perceived  as  being 
too  sharp,  leading  to  artificially  large  differences  between  the  height-wind 
cross-correlations  between  two  residuals  separated  by  just  a  few  degrees  of 
latitude.  This  sudden  change  was  reduced  by  proposing  the  use  of  Bergman's 
expression  (Bq.  3.20)  poleward  of  25*,  then  its  value  at  25*  times  the 
Dey-Morone  expression  (Bq.  BIO)  value  between  10*  and  25*,  with  «  0 
equatorward  of  10*.  A  curve  of  G^  is  shown  in  Fig.  2  as  a  function  of 
latitude,  along  with  the  curve  of  (2kjj)1/2g/fi  for  *  2  x  10-12 
m-1  down  to  25*.  The  latter  curve  is  shown  to  illustrate  how  the  ratio  of 
4  to-  4  ,  used  to  convert  normalized  weights  to  actual  weights  in 

U(V  « 

heights  correcting  winds,  changes  with  latitude.  Between  0*  and  25*  this 
curve  varies  with  pressure  level  because  of  the  tropical  blending  to  form 
4  ;  values  of  4  /4  for  500  mb  are  plotted  as  a  dashed  curve 

Uj  V  U  |  V  4* 

on  Fig.  2.  Using  these  values  of  G^ ,  it  is  seen  that  the  relationships 
between  mass  and  wind  are  primarily  geostrophic  poleward  of  25*,  loosely 
geostrophic  between  10*  and  25*,  and  non-geostrophic  (univariate) 
equatorward  of  10*. 

The  correlations  of  each  selected  residual  with  itself  and  with  all 

other  selected  residuals  are  computed  and  stored  in  a  matrix.  Using 

kl  kl 

Bergman's  symbols,  the  matrix  elements  are  given  by  p^  +  where 

p  is  the  correlation  of  variables  k  and  l's  residuals  at  i  and  j,  and  n  is 
the  correlation  of  observation  errors  of  the  two  residuals.  For  observation 
error  correlation  functions,  Bergman's  relations  are  used  throughout. 
Observational  errors  of  the  same  variable  of  different  instruments  are 
assumed  uncorrelated,  as  are  observational  errors  of  different  sensors  within 

kl 

the  same  instrument  package  (n^  =*  0  for  k  4  1).  The  observational  error 
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Fig.  2  coefficient  of  Geostrophy  and  Ratio  of  Wind  to  Height  Forecast  Error  Standard  Deviation 


kk  kk  kk 

correlation  takes  the  general  form  n  *  y ^ ,  where  y  and  6  represent  the 
horizontal  and  vertical  observational  error  correlations  respectively.  In 
general,  they  take  the  form  of  their  counterparts  y  and  v,  that  is. 


kk  ..  „2  . 

Yij  "  exp<\Sij) 


(9) 


1 

1  +  kp  ln^pj/pj) 


(10) 


Type  1  observations  are  assumed  uncorrelated  between  instruments  in  the 

horizontal,  but  have  a  horizontal  correlation  of  unity  for  the  same 

observation  and  a  vertical  correlation  for  the  same  observation  given 

by  Eq.  10  using  k  =8.3  for  heights  and  5.2  for  winds  (Morone,  personal 
Po 

communication).  Though  winds  measured  by  aircraft  (Type  2)  have  several 
horizontally  displaced  observations  taken  by  the  same  instrument,  no  measure 
of  the  horizontal  spatial  correlation  was  known  at  the  time  of  this 
research,  so  they  are  assumed  to  be  uncorrelated.  Results  from  a  recent  MMC 
report  (Morone,  1986)  on  the  observational  error  correlations  for  aircraft 
winds  could  be  incorporated  into  ASAPOI  in  the  future  to  correct  this 
shortcoming.  Satellite  height  observational  errors  are  assumed  to  be 

-11  -2 

horizontally  correlated  using  Eq.  9  and  k^  =  1.13  x  10  m  and  vertically 

o 

correlated  using  the  values  in  Table  3  (Morone,  personal  communication) 

where  these  values  are  vertically  interpolated  between  values  of  X. 

Correlations  for  satellite  cloud  track  wind  observation  errors  were  unknown, 

so  they  are  assumed  uncorrelated  pending  available  information  on  this  topic. 

kl 

The  right-hand  side  vectors,  made  up  of  the  correlations  for  analysis 
variable  1  =  Z,  u,  and  v  in  turn  with  all  individual  residuals , are  used  with 
the  matrix  to  generate  the  solution  vector  elements  a^separately  for 
1  *  Z,  u,  v.  The  Cholesky  method  (Stobie,  1984)  is  used  to  solve  the  three 
systems  of  linear  equations.  If  a  solution  is  not  found  for  any  reason,  or  if 
one  of  the  normalized  weights  a^  from  any  of  the  three  matrix  problems  is  >  1.1, 
one  residual  from  the  most  highly  correlated  pairs  of  residuals  is  removed, 
the  matrix  is  reduced  by  one  in  rank  and  file,  and  the  problem  re-solved.  At 
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Table  3.  Vertical  Error  Correlation  Function  Values  for  Satellite  T,  z 
(from  Morone,  personal  communication) 


where  X  =  20  (In  (P^Pj)!  +  1 

3  pressure  level  of  location  i 
p,  *  pressure  level  of  location  j 


this  point,  the  actual  weights  a.. ..  are  computed  from  the  normalized  weights 
a^  for  each  of  1  =  Z,  u,  v  analysis  variables  using 


a 

il 


(11) 


where  1  represents  the  variable  being  corrected  by  variable  type  m.  It 

should  be  noted  again  that  4.  is  the  forecast  error  standard  deviation  value 

lm 

for  the  variable  m  from  the  analysis  level  closest  in  pressure  to  the  observation 
*1 

i,  so  that  7“-  may  not  always  be  unity  for  m  =  1.  Also,  in  accordance 
o . 
im 

with  the  use  of  G. ,  G  in  the  observation-analysis  point  correlation  functions 
(G^  is  used  in  a  wind  residual,  height  analysis  point  correlation,  while  G^  is 
used  in  a  height  residual,  wind  analysis  point  correlation),  a^  is  set  to 
zero  for  1  =  Z  when  the  wind  observation  is  equatorward  of  10*,  and  a^  is 
set  to  zero  for  1  =  u,  v  when  the  analysis  point  is  equatorward  of  10*. 

This  insures  that  the  analysis  is  strictly  univariate  in  this  region.  Next, 
the  normalized  analysis  error  (square  of  the  ratio  of  analysis  error  to 
forecast  error  standard  deviation  at  the  grid  point,  see  Bergman's  Eq.  2.13) 
is  computed  using  the  dot  product  of  the  normalized  weight  vector  (for  those 
elements  that  have  non-zero  valued  of  actual  weight  a^  only)  and  its 
counterpart  right-hand  side  correlation  vector.  Then  the  analysis  error  for 
the  analyzed  variable  is  obtained  by  multiplying  the  square  root  of  the 
normalized  analysis  error  by  the  forecast  error  standard  deviation  for  that 
analysis  variable  at  the  grid  point.  This  is  done  for  Z,  u,  v,  but  only 
c»z  is  carried  into  the  next  analysis.  The  correction  for  each  of  Z,  u, 
v  is  calculated  from  the  dot  product  of  the  vector  of  elements  a^  and  the 
vector  of  residuals,  and  its  magnitude  is  checked  against  limits  on  the  size 
of  the  corrections  (250  m  for  Z,  25  m  s-1  for  u,  v).  If  a  computed 
correction  exceeds  its  limit,  the  correction  is  set  to  zero  and  the 
normalized  analysis  error  to  unity. 

Since  temperature  on  the  o  layers  is  the  mass  variable  in  the 
prognostic  model,  the  height  corrections  on  the  layers  have  to  be  converted 
to  temperature  corrections  on  the  layers.  First,  the  01  analysis  is  used  to 
compute  a  height  correction  for  a  a  layer  value  for  a  layer  below 
o  =  1.0  whose  Ao  is  equal  to  that  of  the  lowest  model  layer.  This  is  done 


exactly  the  same  way  it  is  done  for  the  computation  of  height  corrections  in 
all  other  layers,  and  is  done  to  avoid  extrapolation  for  the  lowest  AZ+AT 
conversion.  Next,  the  layer  height  corrections  AZ  are  converted  to  mean 
temperature  corrections  between  the  ML+1  (ML  is  the  number  of  model  layers) 
layer  positions  of  height  correction  using 

"k  --  -‘<AZic+i  -  AV/CR  ln«WV]  •  (12) 

The  mean  temperature  corrections  are  assigned  to  the  "layer-layer”  pressures 
at  the  ML  intermediate  points  between  the  ML+1  a  layer  pressures.  The  a 
values  for  these  intermediate  layer-layer  positions  are  defined  such  that 
the  natural  logarithm  of  the  layer-layer  a  is  the  arithmetic  average  of 
the  natural  logarithm  of  the  a  values  of  the  surrounding  layers.  The 
layer-layer  mean  temperature  corrections  are  then  converted  to  a  layer 
temperature  corrections  using  the  Flattery  algorithm  (from  NMC  global 
spectral  model  preprocessing  code)  described  in  Appendix  A.  This  algorithm 
has  been  used  to  convert  interface  temperatures  to  layer  temperatures,  and 
is  presented  this  way  in  Appendix  A.  Application  of  this  technique  to  the 
problem  of  converting  layer-layer  temperature  corrections  to  layer 
temperature  corrections  is  carried  out  analogously,  and  in  our  case  K  =* 

ML+1,  because  of  the  use  of  a  subsurface  layer.  By  calculating  a  correction 
for  the  subsurface  layer  in  the  upper  air  01  analysis,  the  residuals  are 
extrapolated  to  the  pressure  corresponding  to  using  the  vertical 
structure  functions  in  the  same  way  that  extrapolations  may  occur  at  or 
layers  above  the  highest  nearby  residuals.  A  similar  01  extrapolation  to  a 
layer  above  the  highest  a  layer  was  not  done  due  to  the  sparsity  of  data 
at  those  altitudes.  Uppermost  layer  extrapolation  by  the  Flattery  algorithm 
often  produces  unrealistic  temperature  corrections,  so  for  these  two  reasons 
no  temperature  correction  is  made  in  the  topmost  layer. 

Having  earlier  converted  the  wind  residuals  to  Cartesian  form  for 
analysis  points  poleward  of  70*.  the  computed  wind  component  corrections 
based  on  these  residuals  are  converted  back  to  $,  X  form  before  being 
added  to  the  first  guess  values  to  form  the  analyzed  values.  Temperature 
corrections  from  the  Flattery  routine  are  added  to  first  guess  temperatures 
at  the  analysis  points,  and  analyzed  T,  u,  v  at  all  a  layers  are  written 
on  an  output  file  for  the  analysis  point.  This  is  done  for  all  points  on 
the  61  x  62  analysis  grid. 


c. 


The  Humidity  Analysis  Program  (MASAP) 


Analysis  of  relative  humidity  is  univariate  and  is  performed  in  the 
lowest  MLRH  a  layers.  The  integer  MLRH  is  the  a  layer  index  for  the 
lowest  <j  layer  whose  nominal  pressure  (i.e.,  assuming  p^  =  1000  mb)  is 
less  than  300  mb.  For  the  a  structure  in  this  study,  MLRH  =  7.  MASAP  was 
designed  to  allow  three  options  for  processing  of  humidity  information: 

(1)  regular  01  using  measurements  of  humidity  from  rawinsondes  and  satellite 
soundings,  (2)  regular  01  using  rawinsonde  and  satellite  humidity 
measurements,  supplemented  by  relative  humidity  inferred  from  surface 
weather  observations  (SWO)  and  from  AFGWC  cloud  analyses  (3DN),  and 
(3)  replace  first  guess  relative  humidity  with  3DM  relative  humidity  at 
analysis  points  where  the  cloud  analysis  showed  cloudiness,  and  regular  01 
using  rawinsondes  and  satellite  soundings  elsewhere'.  The  unchecked 
rawinsonde  and  satellite  sounding  relative  humidity  residuals  were  produced 
in  ASAPOI.  For  option  2,  MASAP  performs  similar  processing  to  obtain  SWO 
and  3DN  residuals  in  turn.  All  four  types  are  ordered  from  south  to  north 
with  each  observation  containing  as  many  as  MLRH  residuals. 

After  the  full  set  of  residuals  containing  all  four  types  is  formed, 

the  gross  and  buddy  checks  are  performed.  Forecast  error  standard  deviation 

values  used  for  these  checks  (fixed  in  time)  and  to  begin  an  assimilation 

sequence  were  derived  from  NMC  statistics  for  March  1983  (Morone,  personal 

communication)  and  are  given  in  Table  4.  As  in  the  flagging  procedure  in 

the  ASAPOI  code,  all  observations  are  distributed  into  boxes  of  10“  latitude 

and  of  approximately  equal  area.  Forecast  error  standard  deviation  (FESD) 

values  are  interpolated  to  pressures  of  residuals,  and  residuals  exceeding 

three  times  the  FESD  are  eliminated  in  the  gross  check.  The  buddy  check 

ZZ 

criterion  is  exactly  as  it  is  in  ASAPOI,  where  p. .  is  used  as  the  correlation 

ij 

between  two  residuals,  and  the  lower  quality  observation  of  a  pair  that  fails 

ZZ 

it  receives  a  toss  flag.  If  the  pair  satisfies  the  criterion  and  p„  >  0.75, 
then  if  both  observations  in  the  pair  are  rawinsonde  humidity  residuals  the 
higher  quality  observation  receives  two  keep  flags,  and  if  both  observations 
are  not  from  rawinsondes  the  higher  quality  observation  receives  one  keep 
flag.  QCI  values  for  the  humidity  buddy  check  are  set  as  follows: 
vertically  checked  rawinsondes  have  QCI  =  10,  while  SWO  and  3DN  humidities 
have  QCI  =20.  In  this  way,  rawinsonde  humidities  cannot  receive  toss  or 


Table  4.  Forecast  Error  Standard  Deviation  for  Relative  Humidity 


NMC  Statistics,  March  1983  (Morone,  personal  communication) 


LATITUDE 

Level  (mb) 

90S-10S 

r  (%) 

10S-10H 

r 

10H-30H 

r 

30H-50H 

r 

50H-90H 

r 

1000 

19 

18 

15 

16 

21 

850 

23 

19 

20 

22 

22 

700 

23 

21 

18 

27 

29 

500 

20 

19 

17 

26 

30 

400 

22 

17 

14 

26 

30 

300 

22 

22 

18 

25 

27 

keep  flags  from  other  types,  nor  can  they  impose  keep  flags  on  other  types. 
Satellite  sounding  residuals  can  impose  toss  flags  on  and  receive  keep  flags 
from  SWO  and  3DM,  while  SWO  and  3DM  can  receive  keep  flags  from  and  iopose 
toss  flags  on  other  inferred  residuals.  After  all  toss  and  keep  flags  are 
imposed,  any  residual  having  four  or  more  keep  flags  is  relieved  of  its  toss 
flags  and  is  assured  of  being  considered  in  the  analysis.  All  residuals  are 
then  subjected  to  the  same  type  of  iterative  toss  rejection  procedure  as 
described  for  ASAPOI.  When  this  is  finished,  observations  with  one  or  more 
residuals  remaining  are  assigned  a  number  of  points  for  each  a  layer 
residual:  rawinsonde,  four  points  each;  satellite  soundings,  two  points 
each;  SWO  and  3DH,  one  point  each. 

In  the  01  analysis  for  each  of  the  61x62  analysis  points,  exactly  the 
same  procedure  is  used  to  determine  the  MSMDS  (presently  8)  observations  to 
be  used  in  the  analysis  of  that  point  as  in  ASAPOI.  However,  the  data  types 
in  HASAP  are  considered  in  three  groups:  rawinsondes,  satellite  soundings, 
and  combined  SWO  and  3 DM.  Thus  all  satellite  sounding  humidities  are  only 

ZZ 

considered  if  less  than  NSMDS  rawinsonde  observations  had  >0.1,  and  all 

SWO  and  3 DM  are  considered  only  if  fewer  than  MSMDS  rawinsonde  and  satellite 

ZZ  ZZ 

observations  have  >0.1.  Within  each  type,  *  PTS  is  used  to 
determine  order  of  selection  when  more  than  MSMDS  are  found. 

In  preparation  for  the  computation  of  the  corrections  in  the  univariate 
01,  the  forecast  error  standard  deviation  values  to  be  used  in  the  present 
analysis  are  prepared.  For  the  first  analysis  in  a  data  assimilation 
sequence,  the  values  in  Table  4  interpolated  to  the  analysis  point  latitudes 
and  a  layer  pressures  are  used.  For  subsequent  analyses  in  that 
assimilation  sequence,  the  analysis  errors  from  the  previous  analysis  at 
each  of  MLRH  a  layers  are  read  in,  checked  to  be  >  one  times  the  forecast 
error  growth  amount  but  <  0.40,  and  augmented  by  a  6h  prediction  error 
growth  amount  to  form  the  FESD  values  for  the  analysis  point.  These  growth 
rates  were  estimated  from  early  GSM  runs  and  were  set  at  .08,  .08,  .08,  .08, 
.09,  .10,  .12  respectively  for  the  seven  lowest  NMC  a  layers  from  bottom 
to  top.  At  this  point,  the  relative  humidity  correction  is  set  to  zero  and 
the  analysis  error  is  equated  to  the  just  computed  FESD  for  all  MLRH  layers 
at  the  analysis  point.  These  remain  the  values  for  a  layers  where  no 
subsequent  correction  is  computed  for  lack  of  nearby  observations. 


Processing  for  option  3  takes  place  at  this  point  in  the  code,  in 
advance  of  the  regular  01  analysis.  This  section  is  avoided  if  option  1  or 
2  is  selected.  The  850  mb-300  mb  mandatory  level  values  of  cos[(«/100)*% 
cloud  amount]  from  the  AFGWC  cloud  analysis  at  the  analysis  point  are  read 
in.  If  there  is  at  least  one  of  the  five  mandatory  levels  with  non-zero 
cloud  amount,  the  pressure  level  interfaces  of  the  four  computational  cloud 
layers  (Tibaldi,  1982)  given  by 

p^  a  pA  (the  first  guess  model  surface  pressure) 

P2  =  Px  -50  mb 

p3  =  p2  -(1/3) <p2  -  p5)  (13) 

P4  ■  P3  -(l/3)(p2  -  p5) 

P5  =  300  mb 

are  determined  for  the  analysis  point.  Then  Tibaldi* s  Eq.  4,  which  in  the 
form  used  in  MASAP  is 

rk  *  ^  "  *l  cost(*/100)*%  cloud  amount^]  ,  (14) 

is  used  to  complete  the  inferred  relative  humidity  for  each  of  the  five 
mandatory  levels  indexed  by  k,  where  the  values  of  A^  and  are 
determined  by  which  computational  cloud  layer  L  the  kthpressure  level 
falls  into,  and  the  following  table  from  Tibaldi: 


L 

Cloud  Laver 

% 

4l 

1 

PBL  (?1  -  p2) 

.80 

.20 

2 

LOW  (p2  -  p3) 

.75 

.15 

3 

MIDDLE  (p3  -  p4) 

.60 

.15 

4 

HIGH  (p4  -  p5) 

.55 

.10 

Any  time  any  of  the  levels  falls  below  (greater  in  pressure  than)  p^,  it 
is  assigned  the  PBL  value,  and  any  time  any  level  k  lies  above  (less  in 
pressure  than)  p^,  it  is  assigned  the  HIGH  value.  Then  all  mandatory 
levels  that  receive  a  value  r^  in  this  way  are  used  to  assign  the  same 
value  as  the  a  layer  3DH  by  finding  the  mandatory  level  k  that  is  closest 
in  pressure  to  the  a  layer  pressure.  No  value  of  r  is  assigned  to  a  a 
layer  whose  pressure  is  closest  to  a  mandatory  level  where  r^  was  not 
calculated  due  to  zero  cloud  amount.  Then  for  each  a  layer  that  receives 
an  r  value  in  this  tray,  a  value  of  observation  error  standard  deviation  for 
3DV  (see  Section  IV)  is  interpolated  to  the  a  pressure  and  assigned  as  the 


"analysis  error”  for  that  analysis  point  a  layer  position.  This  direct 
substitution  of  3D1T- inf  erred  humidity  for  first  guess  humidity  constitutes 
the  "analysis”  for  each  cloudy  a  layer  at  the  analysis  point.  Thus,  if 
all  MLRH  a  layers  receive  values  of  r  and  analysis  error  in  this  manner, 
processing  for  that  analysis  point  is  finished.  If  any  o  layer  was  not 
assigned  an  r  by  replacement,  it  is  afforded  an  opportunity  to  receive  a 
Ar  correction  through  the  conventional  01  which  follows.  When  at  least 

zz 

one  observation  has  a  >0.1,  with  a  maximum  of  NSNDS  such  observations, 
corrections  can  be  computed  for  each  of  the  a  layer  positions  at  the 
analysis  point.  The  horizontal  and  vertical  structure  functions  p||T  and 
rr 

vl j  are  identical  in  every  way  to  those  specified  by  Dey  and  Morone  (1985). 
Thus, 


in  which 


rr  rr  rr 

pij  =  Vij  Vij 


(15) 


**5 =  exp  {-vpij/Bg +  "ij11 


(16) 


where  and  ,  the  component  of  distance  between  locations  i  and  j 
parallel  and  normal  respectively  to  the  analyzed  wind  vector  at  the  analysis 
point  and  a  layer,  are  given  by 


1  j — - *±( 'j-'j > 1  (“j' I vg I ) cos  V<tV|VSl)'i"  V 


~ (x, — x . ) ( (+v  / | V  | ) cos  X  +(u  /|V  |) sin  XI) 

o  o  o  &  o  o 


(17) 


M  =  2  f ( X . —x . ) f ( u  /|V  | )cos  X  — (+v  /|V 
ij  1  J  8  8  8  8  8 


)sin  X  ] 
8 


±(y1-y4)t±v_/|V  |)cos  X  +(u  /|V  |)sin  X  J) 


for  !♦  |  >  70*  where  n^,  m^  are  the  map  factors,  x^,  x^,  y^,  y^  are  the 
x,  y  locations  with  respect  to  the  poles  on  a  Cartesian  polar  stereographic 

projection,  u  ,  v  ,  |V  |  are  the  components  and  speed  of  the  analyzed  wind 
8  8  8 

at  the  analysis  point  g,  k-»  2  x  10-12m-2,  B  =  1  +  0.004(u2  +  v2),  and  X.  is 

•*  o  ©  ©  © 

the  longitude  of  the  analysis  point.  The  top  and  bottom  signs  apply  to  the 
northern  and  Southern  Hemispheres  respectively.  Then 


pij  -  <w‘o»«v  V,2l<Y'V,+<v  v<yiv> 

a#  a# 

*ij  =  "  <Vxj)cost(V  ♦j)/21(vg/,vg,)+(V  *j)(ug/,vs,) 

for  |<)>1  <  70*.  Also, 

v"  =  II  +  (kp+  C)ln2(pi/pj)]~1 


(18) 


(19) 


where  k  =  5  and 
P 


-75  3T  ;  3T  <  0  CK  mb1  ) 

dp  dp 

0  ;  dT  >  0 

dp 


(20) 


£  =  -f-  t(Tk+i"Tk)/(pk+r  pk)+<VTk-i>/(V  pk-i)I;  k=2-  MLRH'1* 

Op  4 


(T2-T1)/(p2-  Pl);  k  *  1 


(21) 


-  (VTk-i)/(pk-  pk-i);  k  -  “ 


•  rt* 

where  T^  is  analyzed  a  layer  temperature.  The  correlation  of  each 

residual  of  relative  humidity  at  observation  point  i  with  the  analysis  a 

level  g  are  calculated  using  the  above  equations.  Then  the  quantities 

r»r»  2 

p,  /(l+c.  ),  where  c.is  the  normalized  observation  error,  are  sorted  by  size 
ig  i  i 

and  the  NLRHS  (presently  10)  individual  residuals  with  the  largest  values  are 
selected  for  the  correction  computation.  Here,  observation  error  standard 
deviation  [OESD]  is  interpolated  from  mandatory  levels  to  observation  pressure. 
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but  FBSD  used  is  analysis  level  value.  Then  the  correlation  matrix  elements 
pij  +  nij  cicj  ar®  computed  for  all  pairs  of  the  VLRHS  (or  fewer)  residuals, 
and  the  matrix  equation  involving  these  elements  and  the  right-hand  side 
correlations  p^/(l+c^)  is  solved  for  the  weights  a^.  The  Cholesky  scheme 
is  also  used  here,  with  the  same  provision  made  in  the  case  of  a  singular 
matrix  as  in  ASAPOI.  Then  the  normalised  analysis  error  is  computed  using 
the  weights  and  the  right-hand  side  correlations,  and  the  correction  is 
computed  using  the  weights  and  the  residuals.  Vote  that  since  local 
homogeneity  is  assumed  in  FBSD,  a*  ^  *  a^  since  the  ratio  of  observation 
site  FBSD  to  analysis  point  FBSD  is  unity.  Computed  corrections  must  not 
exceed  0.40,  and  the  resulting  analyzed  value  is  adjusted  (if  necessary)  to 
lie  between  0.1  and  1.  The  analysis  error  passed  to  the  next  analysis  is 
the  product  of  the  FBSD  for  the  analysis  point  and  the  square  root  of  the 
normalized  analysis  error.  Finally,  analyzed  r  is  converted  to  q  using  Bqs. 
1-4  in  reverse  form  using  analyzed  T  and  first  guess  o  layer  pressures. 

rr 

Observation  error  correlations  in  the  humidity  analysis  are  again 
the  product  of  horizontal  (isobaric)  and  vertical  error  correlations. 

First,  it  is  assumed  that  observation  errors  from  two  different  types  of 
observations  (here,  inferred  relative  humidities  are  included  as 
’’observations”)  are  uncorrelated.  It  is  further  assumed  that  all  individual 
rawinsonde  observations  are  uncorrelated.  Vo  attempt  was  made  to  compute 
observation  error  correlations  for  satellite  observations  since  their  OBSD 
values  computed  from  the  FGGB  II-B  data  set  used  were  found  to  be 
excessive.  Values  for  the  horizontal  and  vertical  observation  error 
correlations  for  SWO  and  3DV  are  given  in  Tables  5-8,  having  been  calculated 
following  the  procedure  described  in  Section  IV.  A  SWO  or  3DV  residual  pair 
is  assigned  a  horizontal  correlation  value  according  to  which  50  km 
separation  distance  category  it  falls  into  and  a  vertical  correlation 
according  to  the  closest  mandatory  pressure  level  to  each  residual  in  the 
pair. 

D.  The  Post-Analysis  Program  (POASAP) 

POASAP  converts  the  61x62  analysis  fields  of  T,  u,  v,  q  at  their 
respective  o  layers  into  30  wave  rhomboidal  spectral  coefficients  of 
temperature,  vorticity,  divergence,  and  specific  humidity.  This  is  done 
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Table  5.  Horizontal  Obsarvation  Error  Correlations  for  SUO  Inferred  R.H. 


3 


I 


Distance  (km) 


1000-850 


700-500 


400-300 


5$ 

50-100 

.2039 

.7030 

.9033 

100-150 

.3755 

.5835 

.7196 

\i 

150-200 

.4041 

.5828 

.6305 

200-250 

.3899 

.5363 

.6134 

t 

250-300 

.3620 

.5078 

.5683 

,& 

300-350 

.3035 

.4621 

.5118 

%■ 

350-400 

.3122 

.4515 

.4609 

400-450 

.2912 

.4461 

.4887 

, 

450-500 

.2698 

.3613 

.4395 

500-550 

.2486 

.3596 

.4113 

550-600 

.2334 

.3235 

.3223 

600-650 

.2099 

.2879 

.3558 

$ 

650-700 

.1702 

.2714 

.3600 

4  ‘4 

700-750 

.1859 

.2530 

.3556 

.1 

750-800 

.1675 

.2319 

.2789 

(t. 

l 

800-850 

.1380 

.2053 

.2499 

850-900 

.1421 

.2451 

.3172 

A 

900-950 

.1139 

.1764 

.2583 

•»n 

950-1000 

.1217 

.1984 

.2751 

S* 

1000-1050 

0 

.1629 

.2422 

| 

1050-1100 

.1209 

.2407 

1100-1150 

.1432 

.2301 

1150-1200 

.1440 

.2143 

1200-1250 

.1357 

.1879 

vl 

1250-1300 

.1284 

.2162 

"l*. 

1300-1350 

.1004 

.1734 

$ 

1350-1400 

.1247 

.2075 

1400-1450 

0 

.1588 

••< 

.«* 

* » 

1450-1500 

.1792 

1500-1550 

.2117 

1550-1600 

.1233 

1600-1650 

.1369 

| 

1650-1700 

.1059 

1700-1750 

.1490 

SI 

*  values  in 

this  category  were 

reset  to  0.9 

for  use  in 

humidity  analysis 


Table  6.  Vertical  Observation  Error  Correlations  for  SWO  Inferred  R.H. 


pi-pj 

rr 

Vij 

100O-8S0 

.2867 

1000-700 

.1435 

1000-500 

.0852 

1000-400 

.0843 

1000-300 

.0950 

850-700 

.4727 

850-500 

.3575 

850-400 

.3621 

850-300 

.3496 

700-500 

.5434 

700-400 

.5157 

700-300 

.5084 

500-400 

.8137 

500-300 

.7767 

400-300 

.8895 
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Table  7.  Horizontal  Observation  Error  Correlations  for  3 DM  Inferred  R.H 


Distance  1 

(Ion) 

350 

700-500 

400-300 

0-50* 

1.0 

1.0 

1.0 

50-100 

.2795 

.4967 

.5832 

100-150 

.2121 

.4715 

.5074 

150-200 

.3004 

.4303 

.4917 

200-250 

.2534 

.3941 

.4416 

250-300 

.2604 

.3414 

.4924 

300-350 

.2282 

.3197 

.4055 

350-400 

.1951 

.2752 

.3822 

400-450 

.1928 

.2732 

.4153 

450-500 

.1701 

.2174 

.3814 

500-550 

.1720 

.1892 

.3167 

550-600 

.1586 

.1801 

.3160 

600-650 

.0880 

.1663 

.2908 

650-700 

.1185 

.1152 

.2749 

700-750 

.1167 

.1428 

.3191 

750-800 

0 

.1138 

.2601 

800-850 

0 

.2604 

850-900 

.3270 

900-950 

.2471 

950-1000 

.2744 

1000-1050 

.2712 

1050-1100 

.2378 

1100-1150 

.2270 

1150-1200 

.2277 

1200-1250 

.2499 

1250-1300 

.2043 

1300-1350 

.2333 

1350-1400 

.2453 

1400-1450 

.2431 

1450-1500 

.2303 

1500-1550 

(1750-1800) 

.2112  (.1665) 

1550-1600 

(1800-1850) 

.2444  (.2541) 

1600-1650 

(1850-1900) 

.2091  (.2109) 

1650-1700 

(1900-1950) 

.2405  (.2126) 

1700-1750 

(1950-2000) 

.2199  (.1559) 

*  values  in  this  category  were  reset  to  0.9  for  use  in  relative 
humidity  analysis 


Observation  Error  Correlations  for  SWO  Inferred  R.h. 


pi  -pj 

rr 

Wij 

850-700 

.3833 

850-500 

.1609 

850-400 

.1498 

850-300 

.1346 

700-500 

.4638 

700-400 

.3892 

700-300 

.3863 

500-400 

.7131 

500-300 

.6437 

400-300 

.8026 

using  a  Fast  Fourier  transform  on  the  61  grid  values  on  each  of  the  62 
Gaussian  latitudes,  then  a  Gaussian  quadrature  integration  for  each  Fourier 
component  over  the  latitudes  using  the  appropriate  Gaussian  weights.  The 
analysis  errors  a  ,  a  on  the  61x62  grid  are  similarly  transformed 
to  spectral  form.  Finally,  the  61x62  first  guess  values  of  surface  pressure 
are  updated  at  each  analysis  point  by  extrapolating  analyzed  heights  on  the 
upper  interfaces  of  the  lowest  three  a  layers  (from  analyzed  layer 
temperatures)  down  to  the  terrain  height  Z*.  The  method  used  is  described 
in  Appendix  B.  Then  these  updated  surface  pressures  are  converted  to 
spectral  form  (In  p)  and  included  in  the  set  of  spectral  coefficients  passed 
on  to  the  NMI . 

III.  METHODS  FOR  OBTAINING  REPRESENTATIONS  OF  OBSERVATIONS  ON 
SIGMA  LAYERS/LEVELS 

Much  discussion  among  analysts  of  global  observations  has  focused  on 
the  relative  advantages  and  disadvantages  of  performing  the  analysis  on 
model  a  layers  vs.  isobaric  levels.  Certainly,  the  former  method  has  the 
advantage  of  avoiding  the  necessity  of  vertically  interpolating  first  guess 
(model  forecast)  fields  to  isobaric  levels,  and  then  interpolating  the 
computed  corrections  to  the  a  layers.  It  does  require  an  interpolation  of 
observational  data  to  the  model  a  layers  (or  for  Z,  levels).  The 
motivating  factors  for  using  this  approach  were  (1)  to  avoid  the  smoothing 
of  the  first  guess  fields  and  computed  corrections  resulting  from  the 
vertical  interpolation,  and  (2)  to  make  use  of  all  observations  (including 
significant  level  data)  in  a  sounding  to  minimize  the  error  in  interpolating 
observations  to  a  layers/ levels.  The  following  paragraphs  describe  the 
steps  taken  in  performing  these  interpolations. 

A.  Type  1  Observations  (Rawinsondes ,  Pibals,  Dropsondes,  TWOS  Navaids) 

Once  the  values  for  p,  Z,  T,  DPD,  |v|,  and  0  (wind  direction)  are 
read  in  for  FGGE  II-B  quality  levels  0  or  1  (unchecked,  or  found  valid 
vertical  consistency  check)  for  all  levels  in  a  sounding,  the  observation  is 
checked  to  be  within  three  hours  of  the  analysis  time  and  is  assigned  a 
buddy  check  box.  If  the  observation  is  found  to  lie  within  0.01  degrees 
latitude  and  longitude  of  a  previously  accepted  observation,  the  observation 
with  the  higher  value  of  data  source  index  (DSI)  is  rejected.  DSI  values 
for  various  instruments  in  Type  1  are:  11,  rawinsonde;  12,  pibal;  14,  TWOS 
(tropical  wind  observing  ship)  Navaid;  15,  aircraft  dropwinsonde .  Dew  point 
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depression  (DPD)  given  in  the  observation  data  is  converted  to  specific 
humidity  using  the  expressions 


TD  =  T-DPD 

e  =  6.11  •  exp(19. 9274-5443. 3619/TD) 
w  =  .622  •  e/(p-e) 
q  =  w/(l+w) . 


Wind  components  u,  v  are  computed  from  | V | ,  6  using 

u  =  -|V|sin  6 


v  =  -|V|cos  6  . 


(22) 


(23) 


After  checking  for  duplicate  levels  or  levels  out  of  order  (according  to 
decreasing  pressure),  the  bilinear  interpolation  of  first  guess  T,  u,  v,  r, 
and  pA  to  the  observation  site  is  performed,  a  layer  (p)  and  a  level 

(p)  pressures  are  computed  and  a  level  heights  are  derived  from  ZA  and  the 
bilinearly  interpolated  I  on  o  layers  using  Bq.  5.  If  more  than  one  level 
of  sounding  data  exists  in  the  observations,  the  algorithm  PTOSIG  is 
invoked.  The  method  of  computation  of  a  level  Z,  and  a  layer  u,  v,  T, 
and  q  is  described  in  Appendix  C.  For  Z,  the  observation  levels  immediately 
surrounding  each  a  level  having  valid  values  of  both  T  and  Z  are 
identified,  and  if  they  are  within  300  mb  of  each  other  they  are  used  to 
compute  the  a  level  Z.  Layer  T  is  then  computed  from  each  pair  of  a 
level  Z’s  except  when  one  or  both  is  missing:  then  layer  T  is  missing.  Then 
the  a  level  Z  values  are  recomputed  by  means  of  Eq.  5  for  comparison  with 
their  first  guess  counterparts  using  ZA  and  these  layer  T  estimates.  If  a 
a  layer  T  is  missing  in  the  computation  of  level  Z,  the  Z  value  at  the  top 
of  that  level  is  declared  missing  and  the  first  guess  Z  at  that  level  is 
used  as  the  new  anchor  in  the  computations  for  Z  above.  The  quality  control 
value  assigned  to  each  level  Z  is  the  higher  (poorer  in  quality)  of  the  two 
values  for  the  original  level  Z's  that  led  to  the  layer  T,  each  of  which  in 
turn  was  the  highest  of  the  four  pieces  of  information  (Z^,  Z 2,  T  , 

Tg)  used  from  the  sounding  to  compute  the  original  level  Z  value.  The  300 
mb  difference  limit  is  also  used  for  observation  levels  surrounding  a 
layers  in  the  vertical  interpolation  of  u,  v,  and  In  q.  Specific  humidity 
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interpolation  in  In  q  vs.  In  p  was  performed  to  avoid  the  bias  introduced  in 
a  q  vs.  In  p  interpolation  (Mitchell,  1985).  Vertically  interpolated  q  and 
T  values  on  a  layers  are  used  to  compute  r  on  o  layers  using  Eqs.  1-4. 
Observation-minus-first  guess  residuals  for  Z  (levels),  and  u,  v,  r  (layers, 
with  just  the  lowest  MLKH  layers  for  r)  are  computed  to  finish  the 
processing. 

B.  Type  2  Observations  (Aircraft  Winds) 

Aircraft  observations  consist  of  single  point  observations  of 
temperature  and  wind  at  a  known  flight  altitude.  In  the  PGGE  II-B  data  set, 
several  wind  observations  may  be  included  in  the  report,  each  with  its  own 
latitude  and  longitude  location.  It  was  assumed  that  the  altitude  applies 
to  all  wind  observations  in  the  report. 

The  aircraft  reports,  although  allowing  space  for  it  in  the  record,  do 
not  contain  a  pressure  value  at  flight  altitude.  Because  pressure  is  the 
vertical  coordinate  in  the  system,  this  makes  it  difficult  to  locate  where 
the  observation  should  be  placed  with  respect  to  the  model  atmosphere.  NMC 
(Morone,  personal  communication)  has  provided  formulas  to  assign  a  pressure 
to  an  aircraft  flight  altitude,  as  follows: 


p(Z)  =  226.3  exp (1.576106  X  10  4  (11000-Z)];  Z  >  11  km 
p(Z)  =  1013.5  ((288-0. 0065Z)/288]  ;  Z  <  11  km. 


(24) 


These  formulas  are  apparently  based  on  standard  atmosphere  formulations.  A 

pressure  is  assigned  to  the  flight  altitude  Z  in  the  following  way.  First, 

the  heights  of  the  model  levels  surrounding  the  observation  Z  are  obtained, 

and  the  one  to  which  Z  is  closest  is  identified  (call  this  Z“  ).  Then 

°k 

pressures  are  obtained  for  both  Z  and  Z“  using  the  same  formula  above  (the 

°k 

larger  of  Z  and  Z*  is  compared  with  11,000  m  and  the  appropriate  equation 
°k 

is  used  for  both).  Then  p(Z“  )  is  compared  with  the  model  level  pressure 

°k 

p‘  (Z‘  ),  and  the  difference  is  applied  as  a  correction  to  p(Z). 
a.  a, 
k  k 

At  this  point,  the  model  layer  pressures  surrounding  p(Z)  are 
identified,  and  the  first  guess  temperatures  and  winds  at  chese  layers  are 
linearly  interpolated  in  In  p  to  p(Z).  The  observation-minus-first  guess 
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residuals  individually  for  temperature  and  wind  components  are  calculated  at 
this  level,  and  these  residuals  are  added  to  the  model  layer  first  guess 
values  whose  pressure  is  closest  to  p(Z).  The  respective  sums  constitute  an 
estimate  of  the  observation  value  in  that  model  layer  for  temperature  and 
winds.  A  value  for  height  at  the  top  of  that  layer  is  calculated  from  first 
guess  height  at  the  bottom  of  that  layer  and  this  temperature  using  Eq.  5. 
However,  the  values  of  a  level  heights  derived  from  aircraft  observations 
are  not  currently  used  in  the  ASAP  analysis.  Any  other  wind  observations 
given  in  the  same  report  are  processed  as  above  using  the  same  flight 
altitude. 

C.  Type  4  Observations  (Satellite  Soundings) 

Satellite  sounding  data  in  the  FGGE  II-B  data  set  come  in  four  different 
forms: 

ITYP=1:  layer  thickness  between  reference  pressure  level  and  standard 
pressure  level 

ITYP=2:  layer  precipitable  water  between  a  reference  pressure  level  and 
a  standard  pressure  level 

ITYP=3 :  layer  mean  temperature  between  two  non-standard  pressure 
surfaces 

ITYP=4:  layer  precipitable  water  between  two  non-standard  pressure 
surfaces . 

Satellite  soundings  over  land  are  not  used  in  the  analysis.  The 
discrimination  between  land  and  ocean  soundings  is  made  by  bilinear 
interpolation  of  the  FGGE  Fixed  Field  2.5*  x  2.5*  terrain  field  to  the 
latitude- longitude  position  of  the  satellite  sounding.  If  the  resulting 
height  value  is  greater  than  0.5  meters,  that  observation  is  excluded  from 
the  analysis.  After  the  +3h  time  check  and  assignment  to  a  buddy  check  box, 
the  observation  location  is  checked  against  previously  accepted  observations 
(accepted  Type  1  and  Type  4  obs  before  the  present  one).  If  it  is  located 
within  0.01  degree  latitude  and  longitude  of  an  accepted  observation,  it  is 
rejected.  The  quality  control  indicators  for  Type  4  data  are  set  according 
to  two  criteria:  quality  control  not  made  (QCI=1),  and  report  found  correct 
during  quality  control  (QCI=0).  Other  levels  of  quality  control  are 
included  in  the  data  set,  but  any  observations  not  included  in  one  of  these 
two  categories  are  rejected.  Bilinear  interpolation  of  first  guess  values 
to  the  observation  sito  occurs  next.  Then  the  report  is  checked  to  see  if 
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it  contains  ITYP=1  or  ITYP=3  for  thickness  information  (not  both) ,  and 
ITYPa2  or  ITYP=4  for  moisture  information  (not  both) .  The  following 
paragraphs  describe  the  processing  of  the  four  types  of  data. 

For  ITYP=1  and  ITYP=3,  the  following  procedure  is  used  to  obtain 
representations  of  observed  geopotential  height  on  a  level  surfaces. 

First,  determine  the  observation  layer  mean  temperature.  For  ITYP=3,  this 
is  given  directly.  For  ITYP=1 ,  the  observation  layer  mean  temperature  is 
given  by 


Tk  =  gAZ/ [R  In  (Pk/Pk+1)l;  *=  1.  Kg-1  (25) 

where  AZ  is  the  (positive)  thickness  in  geopotential  meters  between  the 
lower  (in  altitude)  pressure  surface  p^  and  the  upper  pressure  surface 
Pk+1-  Then  for  both  ITYP=1  and  ITYP=3,  the  next  step  is  to  use  the 

Flattery  algorithm  to  estimate  the  satellite  level  temperatures  T.  at  all  K 

1  s 

layer  interfaces  from  the  layer  temperatures  Tfc  One  of  the  resulting  T^’s 
near  the  middle  of  the  sounding  (usually  400  mb)  is  chosen  as  an  anchor 
(TAITCH)  to  recalcu^-ate  the  °th®r  T^'s.  This  is  done  because  it  was  found 
that  the  errors  introduced  by  the  Flattery  algorithm  are  the  smallest  in  the 
middle  of  a  sounding,  with  larger  biases  at  top  and  bottom.  Then  the  other  level 

*  A 

interface  temperatures  T^  are  computed  from  and  the  Tk  using  the  method 

described  in  Appendix  D.  Once  the  T^  are  known  for  all  k,  they  are  used  to 

interpolate  linearly  in  In  p  to  obtain  the  a  layer  temperatures  T  : 

°k 


T,  -  —  CTl+Tu)+{(T)l-Tu)/[ln(p!l/pu)l}{ln  pg  -JLtln(ptPu)  1 )  (26) 

k  2  k  2 


where  T^,  T^  are  the  derived  satellite  level  temperatures  at  levels  p^,  p^ 

immediately  below  and  above  the  o  layer  pressure  p^  .  Then  Eq.  5  is  used 

to  obtain  the  estimates  of  a  level  Z  from  the  T  using  Z*  [or,  at  top  of 

°k 

any  layer  where  T  may  be  missing,  Z.  (o.  ) J  as  an  anchor, 
o.  f .  g .  k 


When  this  method  was  being  developed  it  was  tested  along  with  the 

straight  Flattery  method  using  a  single  rawinsonde  sounding  so  that 

temperature  could  be  used  as  an  input  and  output  heights  could  be  verified 

against  rawinsonde  height.  This  test  is  described  in  Appendix  D.  The 

method  was  also  used  to  compute  RMS  differences  and  biases  in  height  against 

collocated  rawinsondes  for  the  periods  2/3/79  12Z  -  2/22/79  00Z  and  6/13/79 

12Z  -  7/1/79  06Z  using  FGGE  IX-B  data  for  both.  In  this  data  set,  the 

levels  typically  included  in  the  satellite  layer  mean  temperature  soundings 

were  1000,  850,  700,  500,  400,  300,  200,  100,  70,  50,  30,  10,  5,  2,  0.4  mb. 

Using  the  lowest  available  mandatory  level  rawinsonde  height  as  an  anchor, 

the  satellite  layer  temperatures  themselves  were  used  to  estimate  values  of 

satellite  height  at  all  mandatory  levels  above  the  anchor  level  up  to  50 

mb.  The  RMS  differences  and  bias  between  these  and  mandatory  level 

rawinsonde  heights  (except  1000  mb,  the  usual  anchor  level,  and  250  and  150 

mb)  are  given  in  the  left  side  of  Table  9.  Values  on  the  right  side  result 

from  using  reported  satellite  layer  temperatures  to  get  level  temperatures 

according  to  the  method  in  Appendix  D.  Next,  the  values  are  linearly 

interpolated  to  the  mandatory  layer  pressures  where  the  same  relationship 

between  a  level  and  layer  pressures  is  used  to  define  mandatory  layer 

pressures  in  terms  of  mandatory  level  pressures.  Finally,  Bq.  5  and  the 

anchor  rawinsonde  height  are  used  to  estimate  height  on  upper  mandatory 

levels.  The  difference  between  the  two  sets  is  due  to  the  synthesis  of  the 

layer  temperatures,  which  is  done  in  obtaining  o  level  heights  from 

satellite  layer  temperatures  in  the  analysis.  The  RMS  differences  for  the 

synthesized  results  were  used  to  estimate  observation  error  standard 

deviation  values  for  the  satellite  heights.  According  to 

_  1 

*  “22 

Bergman  (1978),  if  the  RMS  differences  [(F  -  F.)  ]  between  observation  sets 

a  d 

F  and  F.  (with  uncorrelated  observation  errors)  are  known,  and  estimates  of 
a  D 

_  1 
2  2 

0ES0  =  (e#)  for  type  a  are  known,  the  same  for  type  b  can  be  evaluated  from 

(e£)  *  -  V{t(F  -  F  )2J  V  -  [(e*)  *]2  .  (27) 

o  a  o  a 
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Table  9.  Geopotential  Height  RMS  and  Bias  for  All  Collocated 
Satellite  Soundings  and  Rawinsonde  Soundings 


*.i  * 


yi 

Vi* 

•ft 


i& 


Using  Reported  Satellite  Layer  T  Using  Synthesized  Satellite  Layer  T 

Level  RMS  (m)  Bias(m)  #OBS  Level  RMS  (m)  Bias(m)  #OBS 
(mb)  (mb) 


02/03/79  122  -  02/22/79  002 


m 

850 

27.13 

-1.55 

293 

850 

24.99 

-2.91 

399 

700 

33.75 

-  .21 

533 

700 

36.23 

-1.11 

481 

1  500 

48.72 

5.70 

561 

500 

50.73 

5.63 

498 

400 

68.55 

10.30 

568 

400 

71.64 

10.83 

501 

»*,V 

$$ 

300 

95.26 

15.72 

559 

300 

99.66 

16.96 

492 

250 

- 

- 

- 

250  * 

98.42 

25.52 

475 

V  1  V 
<*.« 

'A*' 

200 

108.68 

33.82 

521 

200 

102.15 

31.00 

460 

150 

- 

- 

- 

150 

96.43 

24.80 

455 

100 

87.61 

38.11 

490 

100 

85.56 

34.89 

433 

70 

81.25 

36.70 

411 

70 

79.31 

33.22 

355 

50 

91.19 

36.33 

403 

50 

90.16 

31.34 

351 

j| 

,w 

i 

850 

30.55 

-2.61 

06/13/79 

328 

122  -  07/01/79  062 

850  26.62 

-3.36 

502 

f'W 

•s® 

700 

31.13 

.99 

699 

700 

34.02 

-  .13 

645 

500 

30.95 

.81 

746 

500 

33.39 

.26 

678 

400 

52.97 

1.58 

748 

400 

55.64 

1.07 

679 

.  J 

300 

68.14 

1.65 

726 

300 

71.19 

1.72 

660 

/•; 

250 

- 

- 

- 

250 

70.22 

2.41 

652 

200 

61.83 

7.27 

692 

200 

63.68 

3.94 

632 

»{s 

150 

- 

- 

- 

150 

58.28 

-4.36 

620 

100 

66.47 

6.42 

657 

100 

67.07 

.50 

597 

•ty»: 

70 

69.30 

2.15 

585 

70 

73.51 

-3.03 

490 

I 

50 

82.77 

3.19 

565 

50 

87.80 

-1.08 

466 

m\ 


2  2 

Values  from  Dey  and  Morone  (1985)  of  (e&  )  are  used  for  rawinsonde  heights, 
and  RMS  values  from  the  right  side  of  Table  9  are  used  to  produce  estimates 
of  OESD  for  satellite  synthesized  heights  in  Table  10. 

For  ITYP=2  and  ITYP»4,  the  following  procedure  is  used  to  obtain 
representations  of  observed  specific  humidity  on  a  layers.  First, 
determine  the  layer  precipitable  water  for  each  observation  layer.  For 
ITYP*4,  this  is  given  directly.  For  ITYP=2,  we  use  the  integral  nature  of 
the  definition  of  layer  precipitable  water. 


W 


_i_  J  4  q<*p. 
*po  Pu 


(28) 


where  pQ  is  the  density  of  water  and  p^,  py  represent  the  lower 
and  upper  interface  pressures  for  the  layer  over  which  q  is  integrated  to 
compute  the  layer  precipitable  water.  The  integral  over  the  entire  series 
of  layers  is  equal  to  the  sum  of  the  individual  layer  precipitable  water 
values.  Thus  we  may  assume  that  the  V  for  a  particular  layer  is  the 
difference  between  the  U  between  the  top  of  that  layer  and  the  reference 
level  (1000  mb)  and  the  W  between  the  bottom  of  that  layer  and  the  reference 
level.  Now  knowing  the  U  value  for  each  layer  included  in  the  report,  Eq. 

28  may  be  solved  approximately  for  the  value  of  q  at  the  bottom  of  each 
layer  when  an  anchor  value  at  the  top  of  the  uppermost  layer  is  specified. 

In  the  ASAPOI  module  for  satellite  data  processing,  a  value  for  q(pu)  at 
the  top  of  the  highest  layer  in  a  given  observation  is  obtained  by  linear 
interpolation  in  In  p  on  the  first  guess  q  values  to  p^.  Since  this  work 
was  accomplished  before  it  was  discovered  that  In  q  varies  more  linearly 
with  In  p  than  does  q  (see  Mitchell,  1985),  the  specific  humidity  and  not 
its  natural  logarithm  was  used  as  the  dependent  variable  in  the 
interpolation  relationship.  As  satellite  moisture  was  not  included  in  the 
present  assimilation  experiments,  the  bias  was  not  introduced  into  the 
analysis;  however,  future  processing  of  satellite  data  will  involve  the 
assumption  of  linearity  of  In  q  in  In  p.  To  obtain  the  value  of  q(p^) 
at  the  bottom  of  each  successive  layer,  we  assume  linearity  of  q  in  In  p 


between  p.  and  p  : 

1  u 


Table  10.  Computation  of  Observation  Error  Standard  Deviation 
for  Satellite  Height* 


Level 

1 

<e2)2(MMC) 

1 

_  1 

1 

_  1 

A  A  4^  0* 

l<Fr'Fs)  1  <FBB) 

(e2)2(FEB) 

[<F  -F  >2]2<JU1I) 
r  s 

<e2)2(JU*I) 

1000 

7.0 

(12) 

9.75 

(18) 

16.58 

8S0 

8.0 

24.99 

23.67 

26.62 

25.39 

700 

8.6 

36.23 

35.19 

34.02 

32.92 

500 

12.1 

50.73 

49.27 

33.39 

31.12 

400 

14.9 

71.64 

70.07 

55.64 

53.61 

300 

18.8 

99.66 

97.87 

71.19 

68.66 

250 

25.4 

98.42 

95.09 

70.22 

65.47 

200 

27.7 

102.15 

98.32 

63.68 

57.34 

150 

32.4 

96.43 

90.82 

58.28 

48.44 

100 

39.4 

85.56 

75.95 

67.07 

54.28 

70 

50.3 

79.31 

61.32 

73.51 

53.61 

50 

59.3 

90.16 

67.91 

87.80 

64.75 

*  Here,  subscripts  r  a  radiosonde,  s  a  satellite 


u  ].  (29) 


q(p>  =  q(pl)  +  q(Pu)  *  q(V  ~  q(V  tlnp  -  lnpl  *  lnp 


lnp4  -  lnpu 


If  we  substitute  this  expression  into  the  integral  definition  of  W  above,  we 
obtain  the  following  expression  for  q(P|,)>  the  specific  humidity  at  the 
bottom  of  the  layer: 


q<Pt>  • 


W6P0  -  q(Pu) -  Pu>  +  A] 
1  <Pt  -  Pu>  -  A 


(30) 


where 


A(p 


f 


V  = 


pi  - 


u 


In  p,  -  In  p 
l  u 


2  <Pt 


+  Pu) 


Using  the  layer  precipitable  water  W  (in  10-1m)  and  the  top  and  bottom 

_2 

pressures  of  the  layer  (in  mb),  and  using  g  »  9.8  ms  and  p  *  lg 

_3  ° 

cm  ,  q(p^)  will  have  units  of  grams  of  water  vapor  per  gram  of  air. 

If  in  any  lower  layer  the  observed  W  value  is  not  available,  the  computed  q 
value  for  the  bottom  of  that  layer  is  not  available.  A  first  guess  q  value 
is  found  for  that  pressure  level  using  surrounding  qf  values  and  a 
linear  interpolation  in  In  p  to  provide  a  new  anchor  value  for  the 
computation  of  q  from  W  in  layers  below. 

In  the  ASAPOI  code,  values  of  q  at  the  bottom  of  each  report  layer  are 
linearly  interpolated  to  the  o  layer  pressures.  Then  they  are  converted 
to  r  using  Bqs.  1-4  and  the  satellite  c  layer  temperatures  computed  for 
the  same  sounding.  First  guess  a  level  height  and  a  layer  relative 
humidity  are  subtracted  from  their  computed  satellite  counterparts  to  form 
the  residuals  used  in  the  01  analysis. 

A  collocation  study  between  satellite  soundings  and  rawinsondes  was 
conducted  to  determine  values  of  OESO  for  relative  humidity  values  obtained 
using  the  method  described  above.  FGGE  II-B  observations  for  the  same  two 
periods  used  in  the  satellite-rawinsonde  height  study  were  included  in  the 
study,  with  the  typical  soundings  including  values  of  W  for  the  1000-700, 
700-500,  and  500-300  mb  layers.  The  rawinsonde  value  of  q  at  the  top  of  the 
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uppermost  satellite  observation  reporting  layer  (usually  500-300  mb)  was 
used  to  anchor  the  q  calculation.  In  this  way,  a  satellite  observation 
value  for  q  was  derived  for  500  mb  and  700  mb.  Satellite  sounding 
temperature  values  reported  with  the  layer  precipitable  water  values  were 
used  to  convert  q  to  relative  humidity.  At  the  time  of  this  processing,  the 
reported  satellite  layer  temperatures  were  interpolated  from  mean  layer 
pressures  to  the  500,  700,  and  1000  mb  levels  to  perform  this  conversion. 
Since  no  extrapolation  was  allowed,  no  temperature  value  for  1000  mb  was 
computed,  and  thus  no  value  of  satellite  relative  humidity  was  computed  at 
1000  mb.  RMS  differences  for  700  mb  and  500  mb  for  the  February  period  were 
calculated  to  be  34  percent  and  53  percent  in  relative  humidity 
respectively,  while  the  same  calculations  for  the  June  period  resulted  in 
values  of  22  percent  and  25  percent  respectively.  As  the  February  values 
were  thought  to  be  excessive,  scatter  plots  of  collocated  satellite  r  vs. 
rawinsonde  r  were  constructed  to  determine  if  a  correction  relationship 
between  the  two  could  be  derived  empirically.  Because  of  the  significant 
scatter  observed,  no  attempt  was  made  to  try  to  develop  a  correction 
function  for  satellite  r.  Comparisons  of  rawinsonde  r  converted  to  W  with 
satellite  W  were  just  as  disappointing,  with  RMS  differences  in  W  ranging 
from  50  to  70  percent  of  the  averaged  rawinsonde  values.  It  was  decided  not 
to  use  the  Final  FGGB  II-B  satellite  W  observations  in  the  assimilation 
experiments  described  in  this  report  because  of  their  large,  rather  random 
differences  with  rawinsonde  values.  It  is  recommended  that  the 
satellite-rawinsonde  humidity  collocation  study  be  repeated  at  some  later 
date  using  the  FGGE  II-B  reprocessed  data  set,  because  of  the  encouraging 
results  obtained  in  such  a  study  conducted  recently  at  BCMWF  (Lonnberg, 
private  communication). 

0.  Type  6  Observations  (Single-Level  Satellite  Cloud  Motion  Winds) 

Satellite  temperature  and  cloud  motion  wind  data  are  single  level 
observations  like  the  aircraft  data.  In  this  case,  the  pressure  of  the 
observation  p(Z)  is  reported  instead  of  the  altitude.  Therefore,  the 
surrounding  a  layer  pressures  can  be  found  immediately,  and  as  in  the  case 
of  the  aircraft  data,  first  guess  temperature  and  wind  components  are 
linearly  interpolated  in  In  p  to  p(Z).  Observation-minus-first  guess 
residuals  are  calculated  at  this  level,  and  these  residuals  are  added  to  the 
first  guess  values  whose  pressure  is  closest  to  p(Z)  to  form  the  estimate  of 


the  observation  values  at  this  a  layer  pressure.  First  guess  height  at 
the  bottom  of  the  layer  is  used  with  the  computed  layer  temperature  to 
estimate  height  at  the  top  of  the  layer  hydrostatically;  however,  these 
estimated  heights  were  not  used  in  the  present  assimilation  study.  First 
guess  height  and  wind  component  values  are  subtracted  from  their  estimated 
observation  counterpart  to  form  the  residuals. 

IV.  METHODS  FOR  OBTAINING  REPRESENTATIONS  OF  INFERRED 

HUMIDITIES  ON  SIGMA  LAYERS 

A.  Surface  Weather  Observations 

An  algorithm  described  by  Tibaldi  (1982),  modified  slightly  from  a 
method  first  described  by  Chu  and  Parrish  (1977),  was  used  to  estimate  upper 
air  relative  humidity  from  surface  weather  observations  (SWO) .  The  SWO  is 
used  to  obtain  an  estimate  of  the  average  relative  humidity  in  four  layers 
approximating  the  planetary  boundary  layer  and  the  layers  used  for  the 
classification  of  low,  medium,  and  high  clouds  (see  Eqs.  13,  14).  Once  this 
is  done,  the  humidity  at  any  pressure  level  is  set  equal  to  the  layer 
average  humidity  of  whichever  of  these  four  computational  layers  the 
pressure  level  falls  into.  This  is  true  both  for  the  estimation  of  relative 
humidity  at  a  layers  for  the  analysis  and  for  the  values  at  mandatory 
pressure  levels  in  the  collocation  study  to  be  described  later.  The 
following  paragraphs  describe  how  the  four  layer  average  relative  humidities 
were  deduced  from  the  FGGE  II-B  surface  weather  observations. 

After  checking  each  SWO  to  be  within  three  hours  of  the  analysis  time 

(the  SWO  with  the  time  closest  to  the  analysis  time  is  used),  and  making 

sure  that  it  is  not  collocated  with  a  previously  accepted  humidity 

observation  (rawinsonde  or  satellite  sounding),  the  observation  is  checked 

to  see  that  both  the  pressure/height  entry  and  temperature  entry  in  the 

report  have  valid  values.  The  pressure/height  variable  pZ  can  be  given  in 

the  FGGE  II-B  data  set  in  the  form  of  sea  level  pressure,  station  pressure, 

pressure  at  altitudes  ranging  from  1000-4000  m,  or  geopotential  height  at 

850,  700,  or  500  mb.  All  observations  which  have  pZ  in  any  form  other  than 

station  pressure  are  converted  to  station  pressure  using  the  given  values 

p  .  ,  Z  .  from  the  pZ  variable  and  the  inverse  of  the  NMC  "Shuell"  method 
od  ob 


43 


that  is  normally  used  to  reduce  station  pressure  to  sea  level  pressure.  It 
is  then  this  given  or  computed  value  of  station  pressure  that  forms  the 
basis  for  the  definition  of  the  four  computational  layers  whose  bounding 
pressures  are  computed  from  Eq.  13. 

Por  the  boundary  layer,  60  percent  relative  humidity  is  assigned  as  the 

estimate  from  clouds  (RHPBL)  if  either  there  were  no  low  clouds  (C. ) 

L 

present,  or  low  clouds  had  a  base  height  greater  than  or  equal  to  600m.  If 
low  clouds  were  present,  the  value  of  RHPBL  was  set  using  the  expression 
RHPBL=MpBL-ApB£cos[(*/8)0KTASpBLl  in  which  0KTASpBL  had  one  of  two 
values:  (1)  if  in  the  significant  cloud  record,  significant  cloud  types  6 

through  9  are  present  (cloud  types  Sc,  St,  Cu,  Cb),  the  OKTAS  is  set  to 

roL 

the  significant  cloud  amount  for  that  cloud  type,  or  (2)  if  significant 

cloud  data  are  not  available,  but  amount  of  low  or  middle  clouds  (N.  )  is 

b 

nonzero  and  reported  low  cloud  type  is  identified  and  the  cloud  base  height 
is  less  than  or  equal  to  600m,  the  amount  of  low  or  middle  clouds  in  the 
report  is  used  for  OKTASnnr .  Then  the  simple  average  of  RHPBL,  RHTD 
(relative  humidity  from  observed  temperature  and  dew  point  temperature),  and 
RHWW  (relative  humidity  estimated  from  present  weather  table)  is  used  to 
estimate  relative  humidity  for  the  planetary  boundary  layer.  In  this  case, 
QCI=0  is  set  for  the  PBL.  If  one  of  the  three  is  missing,  the  average  of 
the  other  two  is  used  and  QCI=1. 

A  default  value  of  60  percent  relative  humidity  is  also  used  for  RHLC 

for  the  low  cloud  layer  if  no  clouds  are  present.  If  low  clouds  are 

present,  the  value  of  RHLC  was  set  using  RHLC=MLC-ALCcos[(v/8)OKTASLCl 

where  OKTAS,.  is  set  in  the  same  way  as  for  0KTAS__,  except  the  condition 
LC  PBL  r 

for  cloud  base  height  is  removed  in  circumstance  (2)  above.  Then  the  simple 
average  of  RHLC  and  RHWW  is  used  to  give  the  relative  humidity  estimate  for 
the  low  cloud  layer,  and  QCI=0.  If  RHWW  is  missing,  RHLC  is  used  and  QCI=1. 

A  value  of  45  percent  is  assigned  for  the  estimate  of  relative  humidity  in 
the  middle  layer  if  no  middle  layer  (C^)  clouds  (types  Ac,  As,  Me)  are 
present.  If  such  clouds  are  present,  the  significant  cloud  amount  is  assigned 

to  OKTAS  if  the  middle  layer  significant  cloud  information  is  available, 

nc 

and  if  not  the  amount  of  low  or  middle  clouds  (11.  )  is  used  for  0KTAS„_. 

n  nc 

In  either  case,  the  estimate  of  relative  humidity  in  the  middle  layer  is 
obtained  from  M^-A^cosC (w/8)OKTASM(,] ,  and  QCI=0. 


A  value  of  45  percent  is  assigned  for  the  estimate  of  relative  humidity 
in  the  high  layer  if  no  high  layer  (C  )  clouds  (types  Ci,  Cc,  Cs)  are 

H 

present.  If  such  clouds  are  present,  the  significant  cloud  amount  is 
assigned  to  0KTAS„„  if  the  high  layer  significant  cloud  information  is 
available.  If  not,  and  if  the  total  cloud  amount  is  nonzero  and  there  are 
no  low  or  middle  clouds  (N^)  present,  the  total  cloud  amount  (M)  is 
assigned  to  OKTAS„_.  If  the  latter  is  not  true,  high  layer  relative 

HC 

humidity  is  assigned  as  missing.  The  value  of  OKTAS  is  then  used  in 

HC 

MHC-AHCcos1<’/8)0,CTAShc1  to  «9t^inate  relative  humidity  for  the  high 
layer,  and  QCI=0 . 

Once  the  pressure  level  boundaries  and  the  mean  layer  relative  humidity 
for  the  four  computational  layers  have  been  computed  for  each  SWO,  the  a 
layer  relative  humidities  are  assigned  according  to  which  computational 
layer  each  o  layer  pressure  at  the  observation  site  falls  into.  If  any 
o  layer  pressure  is  <  50  mb  greater  than  p  (station  pressure),  it  is 
assigned  the  lowest  computational  layer  relative  humidity.  Any  o  layer 
pressure  <  50mb  less  than  p^  (300  mb)  is  assigned  the  highest 
computational  layer  relative  humidity.  The  quality  control  index  for  each 
o  layer  relative  humidity  is  also  the  QCI  of  the  computational  layer  that 
the  0  layer  pressure  falls  into. 

A  comparison  of  relative  humidity  inferred  from  SWO's  at  mandatory 
levels  and  collocated  mandatory  level  rawinsonde  relative  humidities  was 
conducted  for  both  February  and  June  periods  using  FGGE  II-B  observations. 
For  each  surface  observation  within  0.1*  latitude  and  longitude  of  a 
rawinsonde,  the  observed  surface  pressure  was  used  to  define  the  boundaries 
of  the  four  layers  for  which  the  mean  relative  humidity  was  to  be 
calculated.  Estimates  of  SWO  inferred  relative  humidity  for  mandatory 
levels  between  1000  and  300  mb  inclusive  were  obtained  by  assigning  the 
value  from  the  computational  layer  that  each  mandatory  level  falls  into. 

RMS  differences  and  biases  computed  over  all  observation  sites  and  all  times 
are  given  in  Table  11.  Using  a  OESO  of  5  percent  relative  humidity  as 
suggested  by  Gustafsson  (1981)  for  rawinsondes,  Eq.  27  was  used  to  estimate 
the  0ESD  values  for  SWO  inferred  relative  humidities  based  on  this  February 
study;  these  are  included  in  Table  13  (below).  A  comparison  of  OESO  for  SWO 
relative  humidities  with  the  initial  FESD  values  for  relative  humidity  in 
Table  4  suggests  that  early  in  an  assimilation  series,  the  SWO  relative 


Table  11.  Relative  Humidity  RMS  and  Bias  for  All  Collocated 
SWO  Inferred  Values  and  Rawinsonde  Soundings 


Level  (mb) 

RMS  (X) 

Bias  (X) 

#QBS 

02/03/79 

12Z 

-  02/22/79 

OOZ 

1000 

12.8 

-3.6 

13554 

850 

21.9 

5.1 

23393 

700 

25.2 

3.4 

15986 

500 

24.1 

6.0 

12030 

400 

23.3 

6.9 

9536 

300 

21.6 

5.6 

8457 

06/13/79 

12Z 

-  07/01/79 

06Z 

1000 

12.7 

-3.1 

11876 

850 

21.0 

5.9 

22975 

700 

24.1 

1.6 

14121 

500 

24.9 

7.2 

10870 

400 

24.6 

10.2 

9517 

300 

22.6 

8.4 

8183 

humidities  would  have  about  the  same  influence  on  the  correction  as  do  the 
first  guess  relative  humidities  of  the  6h  forecast.  However,  as  these 
forecast  errors  grow  in  magnitude  in  conventional  data  sparse  regions,  the 
SWO  humidity  residuals  would  have  a  greater  effect  vis-4-vis  the  first  guess 
in  computing  the  correction.  Given  the  reasonable  size  of  the  bias  of  the 
SWO  humidities  with  respect  to  the  rawinsonde  values,  this  could  retard  a 
drift  of  the  analyzed  relative  humidity  field  toward  "model  climatology"  in 
the  assimilation  sequence  as  would  happen  in  these  otherwise  data  void  areas 
without  these  SWO  inferred  humidities. 

The  value  of  5  percent  for  rawinsonde  relative  humidity 
observation  error  standard  deviation  is  probably  optimistic.  Dey  and  Morone 
(1983)  have  suggested  a  value  of  20  percent,  whereas  ECMWF  (II lari,  1986) 
uses  a  value  of  10  percent  at  low  levels  and  20  percent  at  upper  levels. 

The  use  of  5  percent  results  in  significantly  greater  modification  to  the 
first  guess  in  the  01  analysis  than  would  these  larger  values.  However, 
given  the  relative  lack  of  global  rawinsonde  coverage  and  the  tendency  of 
the  model  to  create  its  own  humidity  distribution  over  time  in  an 
assimilation  sequence,  the  5  percent  value  was  used  in  the  present 
experiments.  Values  closer  to  those  used  at  ECMWF  probably  should  be  used 
in  an  operational  environment.  What  is  important  here  is  the  relative 
magnitude  of  the  OESO  for  rawinsonde  relative  humidities  to  that  for  the 
inferred  relative  humidities.  Since  the  emphasis  in  this  study  was  to  use 
inferred  humidities  only  in  rawinsonde  void  areas,  the  indicated  differences 
in  OBSO's  seemed  appropriate. 

The  next  step  was  to  determine  the  spatial  correlations  between 
observation  errors  of  relative  humidity  for  SWO  inferred  values.  In  the 
case  of  inferred  relative  humidity  values,  the  "observation  error”  would 
really  be  a  measure  of  the  error  of  estimate  of  the  relative  humidity 
inferred  from  the  information  given.  It  will  be  assumed  observation  errors 
of  different  observing  systems  are  spatially  uncorrelated,  and  that 
observational  errors  are  uncorrelated  with  the  true  field  values.  Bergman 
(1978)  showed  that,  if  type  ”b"  observations  of  a  meteorological  parameter 
are  compared  with  collocated  type  "a”  observations  in  which  the  latter  are 
assumed  to  be  spatially  uncorrelated  with  each  other,  the  observational 
error  correlation  between  two  type  ”b”  observations  at  locations  i  and  j  is 
given  by 
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or  the  difference  between  observed  values  at 


where  Ar.  =  r.  -  r . .  , 
i  la  ib 

each  point  in  time  and  location,  and  the  overbar  represents  an  ensemble 
average  over  time.  This  means  that  for  the  above  assumptions,  the 
observational  error  correlations  for  observing  system  b  is  equal  to  the 
correlation  of  differences  in  observed  values  between  observing  system  a  and 
observing  system  b. 

These  correlations  were  computed  for  relative  humidity  values  inferred 
from  surface  weather  observations  using  collocated  rawinsondes  as  the  type 
"a**  observation.  Traditionally,  rawinsonde  relative  humidity  observational 
errors  are  assumed  to  be  spatially  uncorrelated,  and  the  same  assumption  was 
made  in  this  study  to  allow  the  correlation  of  difference  to  define  the 
correlation  of  the  alternative  source  estimation  error.  The  above 
correlation  was  computed  for  each  pair  of  relative  humidity  differences  on 
separate  mandatory  levels  over  the  time  period  2/3/79  12Z  -  2/22/79  00Z  (6h 
intervals).  Then  for  each  mandatory  level,  all  correlations  were  grouped  by 
station  pair  separation  distance  categories  (category  1:  0-49.99  km, 

category  2:  50-99.99  km,  ....  category  40:  1950-1999.99  km)  and  were 

averaged  within  each  category,  each  weighted  in  the  average  by  the  number  of 
time  periods  included  in  each  correlation.  Results  of  the  computation  of 
isobaric  correlations  showed  that  the  observation  error  correlations  fell 
into  three  categories:  observations  at  1000  and  850  mb  were  appreciably 
less  correlated  than  those  at  700  mb  and  500  mb,  and  even  less  well 
correlated  than  those  for  400  and  300  mb.  Thus,  the  correlation  values  in 
these  three  categories  were  vertically  averaged  (weighted  by  the  number  of 
pairs  in  all  correlations  over  all  times  at  the  respective  levels).  These 
averages  were  presented  previously  in  Table  5'.  Correlation  values  below  0.1 
are  neglected  in  the  table,  so  that  in  the  case  of  SWO's  distance  categories 
are  included  only  out  to  1750  km. 

Besides  the  computation  of  lateral  correlations,  vertical  correlation 
observational  errors  were  computed  for  SWO  inferred  relative  humidity  at  each 
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station  for  all  categories  of  mandatory  pressure  level  pairs.  These  were 
then  averaged  (each  correlation  in  the  average  weighted  by  the  number  of 
time  periods  included  in  the  correlation)  over  all  stations  to  arrive  at  an 
averaged  correlation  for  each  level  pair  representing  all  stations.  These 
averages  were  presented  previously  in  Table  6. 

B.  3DNEPH  Analysis  of  Cloud  Amount 

A  global  analysis  of  percentage  cloud  amounts,  generated  by  the  Air 
Force  Global  Weather  Central  (AFGWC)  3DNEPH  analysis  model  and  available  in 
magnetic  tape  form,  was  the  source  of  the  cloud  analyses  used  in  this  study 
to  infer  relative  humidity  estimates  on  regular  grids.  Cloud  amounts  are 
available  at  3h  intervals  on  15  levels  in  the  atmosphere  on  one-eighth  mesh 
(~25  nra)  polar  stereographic  grids  (Fye,  1978).  Four  steps  taken  to  get 
estimates  of  relative  humidity  on  mandatory  levels  at  desired  locations 
were:  (1)  determine  cloud  amounts  at  mandatory  levels  (vertical 

compaction),  (2)  horizontally  average  these  to  get  averages  on  one-quarter 
mesh,  (3)  convert  these  cloud  amounts  to  relative  humidity,  and  (4) 
interpolate  bilinearly  to  desired  locations.  In  order  to  perform  a 
collocation  comparison  of  3DNEPH  inferred  relative  humidity  with  rawinsonde 
relative  humidity,  these  four  steps  were  first  applied  toward  creating  files 
of  3DNEPH  relative  humidity  estimates  at  all  rawinsonde  locations  for  the 
period  2/3/79  12Z  -  2/22/79  00Z  at  6h  intervals.  Each  step  is  described  in 
detail  in  the  following  paragraphs. 

In  each  of  the  64x64  boxes  of  points  in  the  512x512  3DNBPH  grid, 
vertical  compaction  was  applied  at  each  eighth-mesh  location.  This  was  done 
by  using  only  the  nine  uppermost  3DNEPH  cloud' layers  in  computing  a  weighted 
average  (by  approximate  pressure  thickness  of  the  layers)  of  layer  pairs  to 
form  estimates  of  cloud  amount  for  the  layers  centered  on  the  mandatory 
levels  850  through  300  mb.  In  particular,  3DNEPH  layers  7  and  8  are 
averaged  and  assigned  to  850  mb,  9  and  10  to  700  mb,  11  and  12  to  500  mb,  13 
to  400  mb,  and  14  and  15  to  300  mb.  In  each  case,  percentage  cloud  amount 
in  the  range  0-100  percent  was  included  in  the  averaging.  An  averaging  was 
done  rather  than  using  the  maximum  since  it  was  important  to  represent  cloud 
amount  through  the  layer  to  estimate  the  layer  relative  humidity  rather  than 
maximum  cloud  cover  at  any  point  in  the  layer.  For  example,  if  3DNEPH 
layers  7  and  8  were  clear  and  overcast  respectively,  a  compacted  cloud 


amount  of  50  percent  cloud  coverage  would  probably  be  more  representative  of 
average  relative  humidity  through  the  layer  centered  on  850  mb  than  would 
the  choice  of  either  the  minimum  (0  percent)  or  the  maximum  (100  percent). 
From  this  point  on  in  the  processing,  mandatory  layers  assigned  zero  cloud 
amount  were  considered  missing  information  for  estimating  relative  humidity. 

The  horizontal  averaging  is  then  performed  on  the  five  compacted 
mandatory  layer  cloud  amount  values  at  each  eighth-mesh  point,  separately  by 
layer.  The  3x3  arrays  of  3DNEPH  points,  each  centered  on  every  other 
eighth-mesh  grid  point  in  both  the  x  and  y  directions,  are  averaged  to 
compute  a  value  of  cloud  amount  for  the  central  grid  point.  This  results  in 
a  quarter-mesh  grid  of  averaged  cloud  amount.  If  all  nine  points  have  cloud 
percent  amounts  between  1  and  100  percent  (no  clear  points),  a  simple 
average  of  all  nine  points  is  completed.  If  some  of  the  grid  points  have 
clear  (no  cloud)  values,  the  eight  outer  grid  points  of  the  3x3  array  are 
divided  into  four  sets  of  triangles,  one  in  each  corner  of  the  3x3  square. 

If  at  least  one  of  the  three  grid  points  in  each  triangle  indicates  the 
presence  of  clouds,  an  average  for  the  center  of  the  3x3  array  is  determined 
by  computing  a  simple  average  of  the  cloud  points.  Unlike  the  vertical 
compaction  in  which  a  layer  average  is  desired,  clear  values  are  not  used  in 
the  horizontal  averaging.  In  the.  case  of  the  horizontal  average  assigned  to 
the  center  point  of  each  3x3  array,  the  desire  is  to  represent  a  cloud 
amount  for  that  point  as  an  average  of  the  cloud  amount  of  the  nine 
contiguous  points,  rather  than  a  cloud  amount  to  represent  an  average  of  all 
of  the  area  included  in  the  3x3  array.  It  was  felt  that  a  cloud  free  3DMEPH 
point  has  no  useful  information  to  contribute  to  estimating  humidity.  If  no 
horizontal  averaging  were  performed,  there  would  be  no  way  to  estimate 
relative  humidity  for  a  cloud  free  grid  point  distinct  from  a  neighboring 
cloud  free  point.  Even  if  a  threshold  humidity  were  assigned  as  is  done  in 
the  Tibaldi  SWO  procedure  at  upper  layers,  all  clear  points  would  receive 
the  same  threshold  value.  In  a  comparison  with  rawinsonde  humidities,  it 
was  felt  that  this  would  increase  the  error  of  estimation  of  relative 
humidity  from  the  3DNEPH  fields  on  the  whole,  although  this  was  not 
confirmed  computationally.  It  is  conceded  that  this  is  just  one  method  of 
horizontally  averaging  the  3DNEPH  cloud  amount  information;  it  is 
recommended  that  others  be  tried  by  comparing  resulting  estimated  humidities 


with  rawinsonde  humidities  to  determine  if  cloud  free  regions  can  supply  any 
useful  information  to  the  estimates.  For  example,  it  may  be  argued  that  at 
the  resolution  of  the  global  analysis,  it  may  be  sufficient  to  average  all 
nine  values  of  any  3x3  array  that  has  at  least  one  cloudy  3DVEPH  point,  and 
think  of  this  as  an  average  over  the  area  of  the  3x3  array  and  assign  the 
corresponding  humidity  inferred  from  this  area  to  any  point  that  lies  within 
it  (i.e.,  the  exactness  of  a  spatial  interpolation  may  be  unwarranted). 

This  would  be  more  in  keeping  with  what  is  done  in  the  vertical  compaction. 

The  next  step  was  to  convert  the  horizontally  averaged  (quarter-mesh 
grid)  and  vertically  compacted  (layers  centered  at  850-300  mb  mandatory 
levels)  cloud  amounts  to  some  humidity  or  functional  cloud  amount  form 
before  bilinearly  interpolating  them  to  the  desired  locations.  This  was 
done  in  three  ways:  (1)  conversion  of  cloud  amount  to  condensation  pressure 
spread  (CPS)  following  the  AFGWC  procedure  (Fye,  1978),  (2)  conversion  of 
cloud  amount  to  square  root  of  the  cloud  fraction  for  use  in  the  ECMWF 
procedure  (Slingo  and  Ritter,  1985),  and  (3)  conversion  of  cloud  amount  to 
cosine  of  (w/100  times  the  percentage  cloud  amount)  for  use  in  the  Tibaldi 
(1981)  procedure.  All  three  of  these  procedures  are  described  in  Appendix 
E.  After  these  conversions,  the  converted  quantities  were  bilinearly 


interpolated  to  desired  locations;  for  the  rawinsonde  collocation  study,  to 
rawinsonde  observation  sites  reporting  at  the  corresponding  synoptic  time. 
Finally,  the  procedures  described  in  Appendix  E  were  used  to  estimate 
relative  humidity  from  the  humidity  or  functional  cloud  amount  forms,  with 
rawinsonde  temperatures  being  used  to  convert  dew  point  depression  to 
relative  humidity  in  the  AFGWC  procedure  in  the  collocation  study. 

RMS  differences  and  biases  for  the  collocated  rawinsonde  relative 
humidities  and  3D1TEPH  inferred  relative  humidities  are  included  for  the 
three  procedures  in  Table  12.  Results  are  presented  for  five  mandatory 
levels  850-300  mb  based  on  all  observation  sites  over  all  6h  time  periods. 

As  in  the  SW0  procedure,  the  Tibaldi  procedure  assigns  the  relative  humidity 
value  to  each  mandatory  level  using  the  coefficients  and  A^  of  Eq.  14 
according  to  which  computational  layer  the  mandatory  level  falls  into. 

Here,  the  observed  surface  pressure  is  used  in  the  computation  of  the 
computational  layer  boundary  pressures.  It  is  clear  from  the  values  in 
Table  12  that  the  Tibaldi  method  results  in  the  smallest  differences  from 
rawinsondes  both  in  terms  of  RMS  and  bias.  Unfortunately,  it  was  not  possible 


Table  12.  Statistics  of  Differences  Between  Relative  Humidity  Inferred 
from  Cloud  Amount  and  That  Measured  with  Radiosondes  (in  X) 


^ /BSSi 

02/03/79 

12Z  - 

02/22/79 

OOZ 

Level 

#0BS 

AFGWC 

KCMWF 

Tibaldi 

(from  Chu-Parrish) 

♦ss®".-- 

tviV 

(mb) 

RMS 

Bias 

RMS 

Bias 

RMS 

Bias 

•sani 

•81 

850 

8251 

20.9 

0.6 

22.2 

9.8 

20.3 

-1.7 

700 

13570 

26.0 

2.2 

29.1 

15.1 

25.3 

-4.5 

m 

500 

10837 

26.3 

-6.8 

32.0 

21.7 

24.1 

-4.7 

400 

4907 

30.0 

-11.6 

34.3 

26.0 

22.7 

-5.4 

300 

7342 

36.8 

-26.5 

37.3 

31.2 

21.3 

-4.7 

to  compute  the  comparable  statistics  for  the  June  period  at  the  time  of  this 
study,  but  it  is  recommended  that  this  be  done  to  verify  that  these  results 
are  not  seasonally  dependent  (because  of  the  greater  number  of  rawinsonde 
observations  in  the  Northern  Hemisphere).  However,  since  the  present 
assimilation  experiments  were  conducted  for  the  February  period,  the  RMS 
values  in  Table  12  from  the  Tibaldi  procedure  were  used  to  compute  the 
observation  error  standard  deviation  values  for  3DN  in  Table  13  using  Eq.  17. 
Horizontal  and  vertical  observation  error  correlations  for  3DN  in  Tables  7  and 
8  were  computed  using  Eq.  31  in  a  manner  directly  analogous  to  that  for  SWO's. 

The  horizontally  averaged  and  vertically  compacted  cloud  amounts, 
converted  to  cosine  of  («/100  times  the  percentage  cloud  amount),  were  also 
bilinearly  interpolated  to  two  regular  grids  for  use  in  the  01  analysis.  For 
moisture  option  2  (described  in  Section  V) ,  values  were  interpolated  to  a 
regular  grid,  3°  latitude  x  120*cos  <t>  points  in  longitude,  to  act  as 
observations  in  the  moisture  01  subroutine.  For  moisture  option  3,  values 
were  interpolated  to  the  61x62  analysis  grid  for  direct  replacement  of  the 
first  guess  relative  humidity  at  cloudy  analysis  points.  In  both  cases,  first 
guess  surface  pressure  (at  the  "observation”  site  and  analysis  point 
respectively)  defined  the  computational  layers  which  dictated  the  choice  of 
M^  and  in  Eq.  14  to  convert  cos[(*/100)*X  cloud  amount]  to  relative 
humidity  for  each  of  the  five  mandatory  level  values  according  to  which 
computational  layer  each  mandatory  level  fell  into.  Then  these  five  relative 
humidities  were  used  to  assign  relative  humidities  to  the  sigma  layer 
pressures  according  to  the  mandatory  level  pressure  lying  closest  to  each 
sigma  layer  pressure. 

V.  DESCRIPTION  OF  DATA  ASSIMILATION  EXPERIMENTS 

In  an  attempt  to  determine  what  effect  each  of  the  three  humidity 
options  has  on  the  humidity  analyses  and  forecasts  initiated  from  them, 
three  five-day  data  assimilation  experiments  were  conducted.  While  it  was 
our  original  intent  to  use  satellite  observations  of  moisture  along  with 
rawinsondes  for  all  01  phases  of  these  experiments,  the  results  of  the 
collocation  study  precluded  their  use,  although  the  software  was  created  to 
include  them.  The  three  humidity  options  are:  (1)  using  rawinsonde 
observations  of  moisture  only,  (2)  using  rawinsonde  moisture,  supplemented 
by  relative  humidity  inferred  from  surface  weather  observations  and  3DNEPH 
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Table  13.  RH  Observation  Error  Standard  Deviation  Values  Used  in 
2/15/79  12Z  -  2/20/79  00Z  Experiments 


Level  (mb)  Rawinsonde  (%)  SWO  Inferred  (%)  3 DM  Inferred  (%) 


*  Assigned  value  not  based  on  collocation  study 


cloud  amounts,  and  (3)  replace  first  guess  relative  humidity  with  3DITBPH 
inferred  relative  humidity  at  cloudy  analysis  points,  and  perform  01  using 
rawinsonde  moisture  at  clear  analysis  points.  The  five  day  period  from 
2/1S/79  00Z  -  2/20/79  00Z  was  chosen  because  of  the  considerable  amount  of 
literature  published  on  comparison  of  analyses  and  forecasts  during  this 
period  [e.g.,  Hollingsworth  et  al . .  (1983)]. 

The  starting  analysis  for  all  three  experiments  was  extracted  from  the 
FGGB  III-A  global  analysis  set  for  2/15/79  at  00Z.  Analysis  values  of  Z,  T, 
u,  v,  on  a  144x73  (2.5*x2.5*  longitude-latitude)  grid  on  mandatory  levels 
1000-50  mb  and  r  on  mandatory  levels  1000-300  mb  were  linearly  interpolated  to 
the  62  Gaussian  latitudes  of  the  analysis  grid,  with  every  other  point  in 
longitude  being  used.  Resulting  72x62  values  of  Z,  T  were  used  in  the  GETPS 
routine  described  in  Appendix  F  to  obtain  corresponding  grid  values  of  pA, 
using  the  corresponding  72x62  FGGB  Fixed  Field  values  of  Z*.  Sigma  layer 
values  of  T,  u,  v  and  q  (lowest  MLRH  layers  for  q)  are  obtained  using  the 
routine  PT0SIG,  where  here  the  data  reside  strictly  on  mandatory  levels.  The 
previously  determined  values  of  pA  are  used  to  define  the  a  layer 
pressures  at  each  point.  Finally,  the  72x62  a  layer  fields  of  T,  u,  v,  q, 
and  the  72x62  values  of  pA,  ZA  are  converted  to  30  wave  rhomboidal 
spectral  forms  of  n,  0,  T,  q,  pA,  ZA  using  the  Fast  Fourier  transform 
and  Gaussian  quadrature  integration  techniques.  Spectral  versions  of  FGGB 
III-A  fields  were  formulated  in  like  manner  at  all  00Z,  12Z  times  in  the 
experimental  period  to  act  as  a  basis  of  comparison  for  the  data  assimilation 
analyses . 

The  Machenauer  non-linear  normal  mode  analysis  was  then  performed  on  the 
spectral  FGGB  III-A  fields  for  2/15/79  00Z,  using  two  iterations  and  four 
vertical  modes.  Then  the  AFGL  global  spectral  model  took  these  initialized 
data  as  input  and  provided  a  6h  forecast  valid  at  06Z.  This  was  the  starting 
first  guess  for  all  three  data  assimilation  experiments.  For  the  first  cycle 
in  all  three  experiments,  the  forecast  error  standard  deviation  values  for  Z 
in  Table  1  and  r  in  Table  4  were  used  in  the  respective  analyses.  In 
subsequent  analyses,  FES0  values  were  obtained  by  augmenting  the  analysis 
error  from  the  previous  analysis  by  the  forecast  error  growth  rate.  Each 
cycle  consists  of  the  sequence  of  four  analysis  steps  (pre-analysis, 
mass-motion  analysis,  humidity  analysis,  post  analysis),  the  normal  mode 


initialization,  and  the  6h  global  forecast.  All  computing  was  performed  on  a 
CDC  CYBER  860  with  367K  octal  words  of  core  memory  and  1720K  octal  words  of 
extended  core  memory.  A  typical  run  of  the  six  programs  in  each  cycle  for 
moisture  option  1  during  00,  12Z  times  took  about  5500  sec  of  central 
processing  time,  of  which  45  percent  consisted  of  the  mass-motion  analysis, 
and  39  percent  was  attributed  to  the  6h  forecast. 

VI.  EXPERIMENTAL  RESULTS 

A.  Mass-Motion  Analyses 

Because  this  series  of  experiments  represents  the  first  test  of  the  ASAP 
data  assimilation  system  in  its  current  complete,  updated  form,  it  was  of 
interest  to  examine  its  ability  to  perform  the  basic  assimilation  of  height 
and  winds.  While  in  the  meteorological  community  there  is  no  definite 
consensus  on  how  to  assess  the  accuracy  of  the  analyses  produced  in  a  data 
assimilation  system,  it  is  agreed  that  the  analysis  should  fit  the 
observations  fairly  well  and  that  this  fit  should  not  change  significantly 
in  time.  Also,  in  a  good  analysis  most  of  the  information  introduced  into 
the  analysis  from  the  observations  should  be  retained  in  the  initialization 
of  the  analyzed  fields.  Finally,  the  longer-term  forecasts  initialized  from 
such  fields  should  verify  with  observations  and  later  analyses  in  a  manner 
that  shows  a  relatively  constant  growth  of  the  forecast  error  with 
simulation  time.  In  the  following  paragraphs,  results  of  the  height-wind 
analyses  will  be  demonstrated  to  show  how  well  the  ASAP  data  assimilation 
system  results  conform  to  these  standards.  In  this  section,  results  from 
Experiment  1  only  will  be  shown  since  the  mass-motion  analysis  was  carried 
out  in  the  same  way  in  all  three  experiments. 

Figs.  3  and  4  show  the  RMS  differences  between  each  of  the  three  fields 
produced  at  each  analysis  time  -  6h  forecast  (GSM),  ASAP  analysis  (ANL) ,  and 

-  -  1  1  2 

the  initialized  analyses  (INT)  -  for  Z  and  V.  Here  V  »  {-  . E.((u-  -  u  )  + 

M  *ai  r  •  o , 

2  1/2  *  * 

(v^  -  vq  )  ]}  for  H  observation  values  uq,  vq  and  field  values  u^,  v^  . 

Corresponding  RMS  differences  for  the  FGGE  III-A  fields,  truncated  to 
rhomboidal  30  waves  and  interpolated  to  a  levels/layers,  are  shown  (F3A) 
for  comparison.  In  fact,  the  F3A  and  ANL  fields  were  used  to  select  the  set 
of  basis  observations  used  in  the  RMS  and  bias  computations.  After 
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computing  observation-minus-field  residuals  for  both  F3A  and  ANL,  both  sets 
of  residuals  were  individually  gross  checked  and  buddy  checked.  Then  only 
those  residuals  that  survived  the  checks  in  both  sets  were  included  as  basis 
observations.  All  other  fields  had  residuals  computed  in  the  same  way,  but 
only  their  residuals  for  the  same  locations  as  the  basis  observations  were 
used  in  the  computations.  The  intent  of  this  procedure  was  to  arrive  at  a 
more  neutral  set  of  observations  against  which  to  compare  the 
field-observation  differences  that  would  not  favor  one  analysis  over  the 
other. 

Immediately  noticeable  from  the  figures  is  the  fact  hat  the 

ANL  values  for  Z  show  a  better  fit  with  observations  than  do  the  F3A 

values  at  all  levels  shown,  while  the  F3A  RMS  differences  for  V  are  smaller 

than  those  for  ANL  at  all  levels  shown  except  o^.  Furthermore,  there  seems 

to  be  no  perceptible  increase  of  RMS  difference  with  time  in  the  ANL  Z  curves, 

yet  with  the  exception  of  the  ANL  curve  all  ANL  curves  for  V  show  somewhat 

of  an  increase  of  RMS  difference  up  to  a  certain  point  and  then  appear  to  level 

off.  An  explanation  for  the  relatively  good  fit  with  observations  for  Z  and 

upper  level  V  and  poorer  fit  for  lower  level  V  lies  in  the  fact  that  these  RMS 

differences  reflect  a  univariate  comparison.  They  are  very  much  dependent 

upon  how  much  influence  Z  observations  have  in  correcting  Z  first  guess 

values,  and  that  u,  v  observations  have  in  correcting  u,  v  first  guess 

values.  Since  both  height  and  wind  observations  are  used  to  correct  both 

height  and  wind  first  guess,  the  fit  of  the  resulting  analysis  with 

observations  will  be  best  in  those  areas  where  observations  of  the  variable 

being  corrected  hold  more  influence  than  their  geostrophic  counterpart. 

Ratios  of  the  normalized  observation  error  for  height  to  those  for  winds, 

shorn  in  Table  2,  play  the  most  important  role  in  determining  the  relative 

influence  of  height  and  wind  observations  surrounding  an  analysis  point. 

Since  normalized  observation  error  affects  the  magnitude  of  the  correlation 
2 

P*  =  p/(l  +  c  )  of  the  observation  with  the  analysis  point, 

observations  with  lower  values  of  c  are  both  more  likely  to  be  chosen  in 

the  analysis  and  will  have  more  weight  in  the  computation  of  the 

correction.  Values  of  c„/c  in  Table  2  show  a  clear  preference 

Z  u,  v 

for  height  observations  over  wind  observations  in  levels  below  200  mb  for 
extratropical  latitudes.  This  agrees  with  the  relative  goodness  of  fit  of 
the  height  analysis  over  the  wind  analysis  at  the  lower  levels.  Notice  that 


at  the  uppermost  levels,  where  rawinsonde  height  errors  are  set  much  higher, 
the  c  /c  ratio  becomes  much  greater  than  unity  and  the  analysis 

2  U  |  V 

is  more  favorable  toward  the  wind  observations.  This  explains  the 
considerable  better  fit  of  the  winds  analysis  to  observations,  which 
occurred  at  o^,  and  a  (~  124,  74,  and  21  mb).  More 
support  for  this  explanation  was  provided  when  the  RMS  fit  with  observations 
was  separated  into  tropical  (J0|<  25*)  and  extratropical  regimes.  In 
the  tropical  latitudes  where  the  analysis  is  much  more  univariate  in  nature, 
the  RMS  differences  between  AML  and  observations  for  V  were  lower  at  all  a 
layers  than  they  were  for  F3A,  whereas  in  the  higher  latitudes  the  F3A  -  ANL 
curves  retained  their  indicated  relationship.  Also,  it  should  be  pointed 
out  that  Bergman  (1979)  used  a  constant  c  /c  ratio  of  1.20  for 

T  U  |  V 

rawinsondes  in  his  analysis,  and  it  was  the  Bergman  analysis  procedure  that 
was  used  to  generate  the  FGGE  III-A  analyses.  This  explains  the  generally 
better  fit  of  the  V  F3A  analyses  and  poorer  fit  of  the  Z  F3A  analyses  as 
compared  with  their  AML  counterparts.  These  results  point  out  the  danger  in 
basing  a  determination  of  quality  of  an  analysis,  especially  one  to  be  used 
to  initialize  a  forecast  model,  solely  on  the  basis  of  goodness  of  fit  with 
observations.  It  is  better  to  use  goodness  of  fit  as  a  trouble  shooter 
criterion  in  analyses  for  initializing  forecast  models,  indicating  problems 
as  they  may  exist  by  showing  significantly  large  RMS  errors  or  large  error 
drifts  in  time  in  a  data  assimilation  mode.  Finally,  a  comparison  of  the 
Z(o^)  RMS  curves  with  the  500  mb  RMS  curves  for  the  NMC  GDAS  (Dey  and 
Morone,  1985)  shown  in  their  Fig.  4  shows  good  agreement  between  the  RMS 

levels,  both  in  the  15-20  m  range  for  the  analysis  and  25-30  m  for  6h 

forecasts.  A  comparison  of  the  V(dg)  RMS  curves  with  the  250  mb  RMS  curves 
in  Fig.  5  of  Dey  and  Morone  shows  that  the  AML,  GSM  curves  appear  to  be 

leveling  off  at  levels  (6.5  -  7.0  ms-1  and  8.5  -  9.0  ms-1  respectively) 

in  fairly  close  agreement  with  the  level  of  the  NMC  curves.  Note  that  since 
the  2/15/79  00Z  starting  analysis  was  the  FGGE  III-A,  the  difference  grows 
from  its  lower  level  throughout  the  five  day  period  to  its  apparent 
equilibrium  level  in  the  lower  o  layers,  while  it  appears  to  decrease  from 
its  FGGE  III-A  level  in  the  upper  a  layers.  This  good  agreement  between 
NMC  and  the  present  AGDAS  results  no  doubt  stems  from  the  close  similarity 
between  the  forecast  models  and  analysis  methods  used  for  both,  and 
especially  from  the  use  of  the  same  observational  errors  for  both  studies. 


In  order  to  verify  that  the  wind  RMS  differences  are  in  fact  leveling  off  at 
the  end  of  the  five  day  period,  the  ASAP  option  1  assimilation  experiment 
should  (but  could  not)  be  carried  out  two  more  days. 

Bias  for  height  and  wind  AGOAS  fields  against  basis  rawinsonde 
observations  was  also  diagnosed  as  a  potential  trouble  indicator;  however, 
only  bias  for  Z  at  four  selected  times  in  the  five  day  period  is 
illustrated  here  (Fig.  5).  The  outstanding  feature  in  the  time  evolution  of 
the  height  bias  is  the  apparently  unabated  growth  of  positive  bias  at  the 
a  =  .050  level  in  the  ANL,  I NT,  and  GSM  fields.  At  all  other  levels,  the 
growth  of  positive  bias  occurs,  but  appears  to  remain  relatively  constant 
after  2/17/79  00Z.  In  all  levels,  the  ANL  appears  to  represent  a  significant 
correction  to  the  GSM  bias,  but  the  ensuing  initialization  tends  to  increase 
the  positive  bias  again.  At  levels  .925  -  .100,  the  fact  that  the  bias  has 
appeared  to  reach  a  steady  value  would  suggest  that  the  analysis  was 
effective  in  compensating  for  the  positive  bias  being  introduced  in  the  6h 
forecast  at  these  levels.  This  is  evidently  not  the  case  for  a  =  .050,  and 
a  possible  explanation  for  this  could  be  the  relatively  large  observation 
error  assigned  at  this  level  (59.3)  resulting  in  relatively  larger 
normalized  observation  errors  and  poorer  fits  to  the  data.  Another  possible 

explanation  is  the  use  of  the  smaller  radial  influence  factor  (1c  =  1.0  x 

-12  -2  " 

10  m  )  resulting  in  larger  correlations  assigned  to  more  distant 

observations  and  thus  a  lessened  fit  with  nearby  observations.  Kistler  and 

Parrish  (1983),  in  a  similar  comparison  of  the  GDAS  with  the  NMC  hemispheric 

nine  layer  model  assimilation,  showed  much  larger  positive  bias  at  100  mb  in 

their  GDAS  in  tropical  regimes,  but  were  more  comparable  in  extratropical 

regimes.  They  attributed  this  difference  to  the  GSM  system's  inability  to 

resolve  the  tropopause  temperature  profile  in  the  tropics  near  100  mb. 

However,  when  the  bias  in  the  ANL,  XNT,  and  GSM  fields  for  2/17/79  00Z  was 

computed  separately  for  latitudes  poleward  and  equatorward  of  25°  latitude, 

A 

the  large  positive  bias  was  evident  in  both  regimes  at  a  =  .050;  in  fact, 
the  positive  bias  was  somewhat  larger  in  the  extratropics.  The  source  of 
this  bias  should  be  studied  further  before  any  really  profitable  analysis 
fields  can  be  produced  at  the  uppermost  levels.  Notice  that  the  FGGE  III-A 
fields  reveal  a  somewhat  more  modest  negative  bias  at  the  a  =  .050  level. 

As  for  the  bias  in  the  wind  components,  results  of  the  five  day 
assimilation  experiment  showed  a  moderate  growth  of  positive  bias  in  the  u 


Observation  Bias  for  a  Level  Heights  (Basis  Observations) 


component  at  a  and  o  and  of  negative  bias  of  v  component  in  wind 

2  1  -1 
at  o  ,  both  of  which  appear  to  level  off  at  about  the  1.0  ms  level 

O 

about  mid-way  through  the  period.  Also  of  interest  was  a  decrease  in  a 

positive  bias  in  u  at  o .  in  the  GSM,  INT  fields  only  from  about  1.8  m 

-1  -1 
s  at  the  beginning  to  1.2  m  s  at  the  end  of  the  period. 

Interestingly,  no  significant  bias  was  apparent  in  the  corresponding  ANL 

fields.  Thus  it  appears  that  the  wind  biases  suggest  no  apparent  trouble 

areas  in  the  wind  data  assimilation. 

Time  evolution  of  the  global  average  of  the  analyzed  fields  was 

examined  as  an  additional  indicator  of  meteorological  realism  of  the 

analyses.  Figs.  6-8  illustrate  the  change  of  global  average  of  temperature, 

kinetic  energy,  and  surface  pressure  through  the  five  day  period.  Global 

averages  of  the  FGGE  III-A  fields  are  shown  for  comparison.  While  no  clear 

trends  are  evident  in  the  temperature  curves  for  a .  and  o, ,  both 
-  4  7 

02  and  curves  show  a  clear  growth  of  the  globally  averaged  a 
layer  temperature  with  time.  At  the  lower  level,  there  appears  to  be  a 
leveling  off  of  the  growth  after  an  increase  of  about  0.8°K  at  2/18/79  00Z; 
however,  this  should  be  confirmed  through  following  analyses.  Mo  such 
indication  of  leveling  off  is  evident  in  the  curve,  and  in  fact  the 
warming  indicated  at  that  layer  is  consistent  with  the  growth  of  positive 
bias  of  Z  at  a  *  .050  (the  level  at  the  top  of  o  layer  11)  as  shown  in 
Fig.  5.  The  global  mean  temperature  grows  so  fast  in  time,  as 
compared  to  the  slight  trend  toward  cooling  at  o^,  the  o^0  layer 
becomes  the  layer  of  minimum  global  mean  temperature  (labelled  "TROP")  at 
2/18/79  00Z.  The  F3A  o^,  curves  maintain  the  relative  global 
average  layer  temperature  relationship  for  these  two  layers.  This  gives 
additional  support  to  the  idea  that  the  anomalous  warming  in  the  data 
assimilation  period  is  confined  to  the  layer  above  a  =  .100.  Fig.  7  does 
not  reveal  any  similar  troubling  features  in  the  evolution  of  the  global 
average  kinetic  energy.  There  appear  to  be  similarities  between  the  trends 
reflected  in  the  F3A  and  ANL  curves.  For  example,  both  oQ  curves  show  a 
growth  of  kinetic  energy  through  2/17/79  12Z,  then  a  relative  leveling  off 
of  the  trends  beyond  that  point.  The  surface  pressure  curves  in  Fig.  8 
reveal  no  untoward  growth  or  decline  through  the  period.  F3A  curve 
indicates  a  small  trend  toward  an  increase  in  the  global  average  surface 
pressure,  while  the  AML  curve  suggests  an  even  smaller  tendency  toward 
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Fig.  6  Global  Average  a  Layer  Temperatures  for  Analysis  Fields 
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decrease,  though  probably  statistically  insignificant.  The  accompanying  a 
layer  temperature  curves  in  Fig.  8  are  shown  because  of  the  fact  that  these 
are  used  with  surface  geopotential  to  obtain  a  level  geopotential  at  the 
lowest  three  a  levels,  from  which  analyzed  surface  pressure  is  diagnosed 
(see  Appendix  B) .  It  is  interesting  that  the  layer  temperature  increases 
seen  in  the  first  half  of  the  assimilation  period  are  not  reflected  in  any 
perceivable  way  in  the  diagnosed  surface  pressure.  Further  study  of  the 
surface  pressure-analyzed  geopotential  relationship  should  be  conducted  to 
better  understand  the  effect  of  change  of  geopotential  on  change  of  surface 
pressure. 

Another  measure  of  the  effectiveness  of  an  analysis  method  in  imparting 
useful  observational  information  to  the  correction  of  the  forecast  first 
guess  is  the  amount  of  observational  information  retained  in  the  analysis 
after  it  has  been  initialized.  In  removing  the  fast  gravity  mode  tendencies 
in  the  analysis  field,  the  initialization  procedure  produces  a  field  which 
contains  only  the  slow  modes  (lower  frequency  oscillations),  which  are  the 
only  ones  useful  to  the  model.  It  is  then  of  interest  to  determine  how  much 
observational  information  imparted  by  the  analysis  consists  of  slow 
modes — that  is,  how  much  is  left  behind  after  the  analysis.  The  RMS 
difference  between  the  analyzed  and  initialized  fields  is  a  measure  of  the 
fast  mode  information  in  the  analysis  that  was  removed  due  to  the 
initialization.  Similarly,  the  RMS  difference  between  the  forecast  and 
initialized  forecast  is  a  measure  of  the  fast  mode  information  introduced  by 
the  forecast  model.  The  difference  between  the  second  difference  (fast  mode 
in  forecast)  and  the  first  difference  (fast  mode  in  analysis)  would  then  be 
the  fast  mode  information  resulting  from  the  update  to  the  forecast  field  by 
the  analysis.  Then  if  the  initialized  analysis  contains  only  slow  modes  in 
the  analysis  field,  and  the  initialized  forecast  contains  only  slow  modes  in 
the  forecast  field,  the  RMS  difference  between  the  two  would  be  a  measure  of 
the  slow  modes  due  to  the  analysis.  Values  for  the  analysis  fast  mode/slow 
mode  ratio  based  on  RMS  differences  averaged  over  the  ten  time  periods  along 
with  the  RMS  differences  between  the  fields  themselves  are  given  in  Tables 
14  and  15.  Lower  values  of  this  ratio  indicate  a  greater  effectiveness  of 
the  analysis  method  in  imparting  the  useful  slow  mode  information  into  the 
corrected  first  guess.  Results  show  a  clear  increase  of  this  ratio  with 
altitude  for  heights  but  a  relatively  steady  value  with  altitude  for  winds. 


Cl:  I NT  GSM 


Unfortunately,  it  was  not  possible  to  obtain  this  ratio  for  FGGB  III-A  as  a 
comparison  as  all  four  assimilation  fields  are  necessary  for  this 
computation.  It  is  apparent  that  at  upper  levels  for  heights,  the 
initialization  makes  a  significantly  larger  change  in  the  analysis  (over  90 
percent  of  the  G-I  difference  at  than  at  lower  levels  (approximately 

A 

10  percent  of  the  G-I  difference  at  02 )•  Apparently  the  analysis  provides 
a  considerably  lesser  amount  of  useful  observational  height  information  to 
the  forecast  correction  at  upper  levels  than  at  lower  levels. 

To  the  weather  forecaster,  the  best  measure  of  an  analysis  is  the 
accuracy  of  the  forecast  produced  from  it.  The  closeness  of  fit  of  the 
analysis  or  the  initial  conditions  is  of  secondary  importance  to  the  fit  of 
the  forecast  field  to  observations  at  the  valid  forecast  time.  As  a  test  of 
the  conformity  of  the  ASAP  analyses  to  this  standard,  a  48h  forecast  using 
the  GSM  was  conducted  using  the  initialized  FGGB  III-A  and  ASAP  analyses  at 
2/17/79  12Z,  exactly  half-way  through  the  five  day  assimilation 
experiments.  Resulting  spectral  forecast  fields  at  12h  intervals  were 
evaluated  at  1*  x  1*  grid  intervals,  and  bilinearly  interpolated  to  the 
locations  of  the  basis  observations  for  each  verification  time.  Meanwhile, 
the  observations  were  interpolated  to  the  0  level/ layer  pressures  of  the 
forecast  fields  and  the  observation-minus-forecast  field  residuals  were 
computed.  The  basis  observation  subset  of  these  residuals  was  used  to 
construct  the  RMS  and  bias  curves  shown  in  Figs.  9-13.  The  RMS  and  bias 
curves  for  the  2/17/79  12Z  analyses  are  included  in  the  first  plot  as  a 
reference.  The  RMS  curves  reveal  an  interesting  relationship  between  the 
importance  of  the  fit  of  the  initial  height  fields  to  observations  and  the 
initial  fit  of  winds  in  the  fit  of  the  forecast  field.  The  better  fit  of 
the  ANL  (ASAP)  analysis  initially  is  soon  lost  in  the  forecast,  especially 
at  the  upper  levels.  The  gap  between  the  height  forecast  RMS  for  FGGB  III-A 
and  ASAP  widens  gradually  with  time  at  those  levels,  to  become  about  12  m  at 
0  3  .100  for  48h.  However,  by  this  time  differences  below  about  500  mb  are 
slight.  By  contrast,  the  gaps  between  the  V  RMS  curves  remain  essentially 
constant  in  time,  widening  only  slightly  from  12h  to  48h.  It  appears  that 
the  initial  fit  played  a  more  significant  role  in  the  forecast  verification 
for  the  winds  than  for  the  heights.  All  eight  sets  of  curves  move 
constantly  toward  higher  RMS  values  with  time,  but  the  wind  curves 
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Differences  on  a  Levels  between  Forecast  Height  Fields  and  Basis  Observations 


Differences  on  a  Layers  between  Forecast  Vector  Wind  Fields  and  Basis  Observations 


Observation  Bias  for  a  Level  Heights  (Basis  Observations) 
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Fig.  12  Forecast  Field  -  Observation  Bias  for  a  Layer  Zonal  Component  of  Wind  (Basis  Observations) 
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essentially  retain  their  relationship  whereas  the  gap  between  the  height 
curves  widens.  Caution  must  be  taken  in  attempting  to  draw  too  many 
conclusions  from  a  single  forecast  experiment.  A  series  of  forecast 
comparisons  should  be  conducted,  especially  for  later  initial  times  when  the 
analyzed  wind  fields  have  attained  their  steady  values  of  RMS  fit  with 
observations  (say,  not  before  2/20/79,  00Z) . 

The  height  bias  curves  in  Fig.  11  reveal  two  major  trends  in  the 

forecast  for  FGGE  III-A  and  ASAP  initial  fields:  (1)  a  positive  bias 

increasing  steadily  in  time  at  all  levels,  and  (2)  a  faster  growth  of 

positive  bias  in  the  ASAP  initialized  forecast  at  level  a  =  .150  and  above, 

and  a  faster  growth  of  the  positive  bias  in  the  FGGE  III-A  initialized 

forecast  at  level  a  *  .375  and  below  after  12h  and  at  level  a  =  .250  and 

below  after  24h.  Particularly  striking  is  the  widening  of  the  gap  between 

the  two  curves  at  a  =  .100  and  .050.  This  is  primarily  due  to  the  original 

significantly  large  positive  biases  at  these  levels  in  the  initial  analysis 

for  ASAP,  and  the  large  negative  bias  in  the  FGGE  III-A  initial  analysis  at 

a  =  .050.  The  forecast  model  appears  to  take  these  large  initial 

differences  and  increase  them  even  more  during  the  course  of  the  forecast. 

This  widening  at  these  uppermost  levels  is  due  primarily  to  the  extreme 

growth  of  the  ASAP  initialized  positive  bias.  A  comparison  of  Figs.  5  and  11 

for  the  AML  (ASAP)  curves  would  indicate  that  much  of  the  growth  of  the 

* 

positive  bias  at  the  a  =  .050  level  in  the  assimilated  analyses  (Fig.  5) 
must  be  due  to  the  characteristic  of  the  GSM  to  create  a  strong  positive  bias 
at  this  level,  as  opposed  to  the  model  used  to  assimilate  the  FGGE  III-A 
analysis  which  left  a  modest,  fairly  constant  negative  bias  at  this  level. 

The  wind  component  bias  curves  (Figs.  12  and  13)  are  similar  to  the 
vector  wind  RMS  curves  (Fig.  10)  in  that  they  retain  their  relationship  with 
each  other  quite  well,  especially  after  12h.  Looking  first  at  the  u 
component  bias,  we  see  in  both  curves  a  tendency  toward  positive  bias  in  the 
lower  layers  (with  ASAP  based  forecast  being  more  positive)  and  toward 
negative  bias  in  the  tropospheric  extratropical  jet  and  levels  immediately 
above  it  (with  FGGE  III-A  based  forecast  being  more  negative).  The 
extremely  high  positive  bias  at  the  uppermost  layer  for  FGGE  III-A  based 


forecast  is  present  in  the  initial  analysis  and  is  continued  in  the  model 
forecast.  It  is  due  to  assigning  the  FGGE  III-A  analysis  50mb  winds  to  the 
a  *  0.21  (~21mb)  level,  which  in  the  case  of  the  zonal  wind  would  lead 
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to  an  overestimate  as  is  evidenced  by  the  positive  bias.  This  direct 
assignment  was  used  in  lieu  of  extrapolation,  because  50mb  is  the  highest 
level  carried  in  the  FGGE  III-A  analysis  fields.  Interestingly,  the  low 
layer  positive  bias,  high  layer  negative  bias  is  also  present  in  the 
meridional  component  forecast  (Fig.  13).  Here,  however,  the  values  of  bias 
between  the  forecasts  are  more  similar,  especially  at  48h.  The  FGGE  III-A 
based  forecast  does  exhibit  a  tendency  toward  larger  negative  bias  just 
below  the  tropospheric  extratropical  jet  level  (in  a  layers  .436  -  .275). 

B.  Humidity  Analyses 

For  the  univariate  relative  humidity  analysis,  fit  of  analyses  to 
observations  becomes  somewhat  more  important  since  only  humidity  data  are 
involved.  In  this  case,  the  RMS  fits  are  very  sensitive  to  the  observation 
error  levels  assigned  to  each  observation  type.  This  is  illustrated  in 
Figs.  14  and  15,  which  show  the  fit  of  the  ASAP  analysis  and  forecast  fields 
to  the  basis  observations  for  moisture  option  1  (rawinsonde  humidity  only 
used  in  the  analysis).  Not  surprisingly,  the  use  of  5  percent  OESO  for 
rawinsonde  humidity  observations  results  in  a  closer  fit  with  those 
observations  than  does  the  use  of  20  percent  as  in  the  F3A  for  all  layers. 

The  ANL  curves  show  lower  RMS  differences  and  biases  and  indicate  no 
apparent  drift  away  from  observations.  The  sawtooth  configuration  of  the 
GSM  bias  plot  for  appears  to  have  a  diurnal  character  to  it,  while 
the  highly  irregular  bias  curves  at  a ?  are  probably  due  primarily  to 
fewer  comparison  observations  for  that  layer. 

Fig.  16  shows  the  distribution  of  checked  rawinsonde  observations 
summed  over  10*  latitude  bands  and  over  the  ten  00Z,  12Z  observation  times 
during  the  five  day  period.  Clearly,  the  first  guess  remains  relatively 
uncorrected  south  of  about  50*  south  latitude,  with  the  maximum  corrections 
in  conventional  01  relative  humidity  analyses  occurring  in  lower 
tropospheric  levels  in  Northern  Hemisphere  midlatitudes.  Fig.  17(a-c)  shows 
the  zonal,  time  average  cross  section  of  (a)  observed  relative  humidity 
using  the  checked  observations  whose  distribution  is  depicted  in  Fig.  16, 

(b)  the  analyzed  (moisture  option  1)  relative  humidity  resulting  from  the  ASAP 
analyses  averaged  over  all  00Z,  12Z  times  during  the  five  day  assimilation, 
and  (c)  the  results  of  subtracting  plot  (a)  from  plot  (b);  i.e.,  the 
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Pig.  16  Total  Humber  (XIO)  of  a  Layer  Rawinsonde  Values  of  Relative 

Humidity  Summed  Over  10*  Latitude  Bands  and  the  Ten  00.12Z  Time 
Periods  of  the  Data  Assimilation  Experiment.  These  Are  the 
Checked  Observations  Making  Op  the  Set  of  Basis  Observations. 
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Fig.  17(a-c)  Zonal,  Tim*  Average  Ralativa  Humidity  for  (•)  Basis 

Observations ,  (b)  ASAP  (Opt.  1)  Analyses,  and  (c)  Difference 
between  (b)  and  (a).  Contour  Interval  for  (a),  (b)  is  5%, 
for  (c),  2%.  Dashed  Contours  Indicate  Hegative  Values. 
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differences  in  the  zonal,  time  averages  of  observed  and  analyzed  a  layer 
relative  humidities.  The  a  layer  pressures  indicated  on  plot  (a)  are  the 
nominal  (pA  =  1000)  layer  pressures,  whereas  those  shown  in  plots  (b)  and 
(c)  are  based  on  the  ASAP  weekly  averaged  surface  pressure.  The  dense 
contour  intervals  on  plots  (a)  and  (c)  are  the  result  of  very  few 
observations  lying  in  the  latitude  band  90*S  -  80*S.  The  difference  plot 
(c)  reveals,  in  a  zonal,  time  average  sense  (ignoring  the  left-most  quarter 
of  the  plot) ,  a  generally  good  fit  of  the  observations  especially  between 
the  848  mb  and  332  mb  levels.  There  is  a  tendency  for  the  analyses  to  be 
too  moist  at  the  lowest  layer  for  most  latitudes,  especially  those  dominated 
by  oceans.  This  is  probably  a  result  of  the  forecast  model  in  that  moisture 
fluxes  from  the  ocean  surfaces  are  not  being  mixed  into  upper  levels  by  the 
rather  simplistic  boundary  layer  parameterization.  The  moist  bias  also 
appears  at  most  latitudes  for  the  highest  (270  mb)  layer,  though  the  reason 
for  this  is  not  as  clear.  However,  with  the  exception  of  under-analyzing 
the  dryness  of  the  descending  branch  of  the  Hadley  circulation  at  about 
15°H,  the  fit  at  other  layers  is  quite  satisfactory.  A  comparison  of  Fig. 
17(c)  with  the  corresponding  difference  cross  section  for  FGGE  III-A 
analyses-minus-observations  (Fig.  18)  reveals  once  again  the  effect  of 
lowering  the  OESO  value  on  the  analysis  fit  to  observations.  Part  of  the 
increased  moist  bias  in  the  FGGE  IXI-A  difference  cross  section  is  due  to 
the  necessity  of  extrapolating  humidity  from  the  highest  analysis  level  (300 
mb)  to  the  a ^  layer  pressure  (usually  less  than  300  mb).  However, 
through  the  lower  layers  a  pronounced  moist  bias  is  evident,  whereas  at  a 
layer  pressures  565,  430,  and  332  mb  a  definite  dry  bias  in  the  analysis  is 
indicated.  These  results  confirm  in  a  zonal,  time  average  sense  what  Fig. 

15  showed  in  a  synoptic,  point-by-point  verification.  The  F3A  biases  are 
not  present  in  the  ASAP  analyses,  probably  largely  due  to  the  lower  error 
assigned  to  the  observations. 

It  is  in  the  data  sparse  regions  as  indicated  in  Fig.  16  that  the  use 
of  alternative  forms  of  moisture  information  would  be  expected  to  have  their 
greatest  effect.  This  can  be  illustrated  in  one  way  by  comparing  the  plots 
of  observation  locations  for  moisture  option  2  (rawinsonde,  SWO,  and  3DV 
relative  humidities  used  in  the  analysis)  as  shown  in  Fig.  19(a-c).  These 
figures  show  the  distribution  of  (a)  unchecked  observations,  (b)  checked 
observations,  and  (c)  observations  actually  used  at  least  once  in  the 
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Fig.  18  Difference  between  Zonal,  Time  Average  FGGK  III-A  Relative  Humidity 
Analyses  and  Zonal,  Time  Average  Relative  Humidity  Basis 
Observations.  Contour  Interval  is  2%,  and  Dashed  Contours 
Indicate  Regative  Values. 
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Pig-  19(a,b)  Location  of  (a)  Unchackad  and  (b)  Chackad  Ralativa  Humidity 
Obaarvations  for  ASAP  (Opt.  2)  Analysis  of  2/15/79,  12Z  at 
o,  (-860  ab).  R  ■  Rawinsonda,  Z  -  Surface  Waathar 
Observation  Inferred,  •  -  3DTOPH  Inferred. 
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Pig*  19c  Location  of  Relative  Humidity  Observations  Actually  Used  in 

2/15/79,  12Z  ASAP  Analysis.  Symbols  Are  Same  As  in  Fig.  17(a,b). 
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analysis  for  2/15/79  12Z  at  a (-860  mb).  Observations  present  in  (a) 
but  missing  from  (b)  were  removed  by  the  buddy  and  gross  checks.  By 
comparing  plots  (b)  and  (c),  one  can  see  the  dominance  of  the  rawinsonde 
observations  being  selected  over  alternative  types,  even  though  the 
rawinsonde  observations  may  be  many  fewer  in  number  and  more  scattered  in  an 
area.  For  example,  over  India  a  predominance  of  SWO  humidities  (and  to  a 
lesser  extent,  3D1I  humidities)  exists  in  the  checked  data,  but  the  analysis 
preferred  to  use  the  fewer  and  more  scattered  rawinsonde  observations  for 
analysis  points  in  that  region.  Alternative  sources  are  all  but  excluded 
from  use  over  rawinsonde  dense  regions.  This  is  due  to  the  much  lower 
observational  error  for  rawinsondes  which  favored  their  selection  and 
resulted  in  higher  weights.  Thus,  major  differences  between  option  1  and 
option  2  humidity  analyses  should  exist  only  in  the  rawinsonde-sparse 
regions,  such  as  0*-30*N  and  60*-90*H.  It  was  thought  that  by  comparing 
analyses  against  zonal,  time  averages  of  observations  that  the  data  scarcity 
in  those  regions  could  be  overcome  to  allow  for  reasonable  comparisons . 

Plots  of  zonal,  time  average  option  2  analyses  and  their  difference  from 
zonal,  time  average  basis  observations  (i.e.,  rawinsondes)  are  given  in  Fig. 
20(a-b).  Indeed,  we  see  the  greatest  differences  between  Figs.  17(c)  and 
20(b)  south  of  30*11,  where  in  the  mid-layers  (848  mb  -  429  mb)  the  prominent 
differences  are  two  locations  of  greater  moist  bias  in  the  Hadley  branch 
(around  15*10  and  centered  at  the  712  mb  layer  at  20*-30*  S.  The  option  2 
analyses  appear  to  be  significantly  more  moist  at  these  layers  and  levels, 
and  because  of  the  conventional  data  sparsity  in  these  regions  the 
differences  must  be  due  to  the  use  of  the  alternative  moisture  sources. 
However,  note  that  a  positive  effect  of  using  the  alternative  moisture 
sources  is  that  the  moist  bias  is  lessened  at  the  lowest  and  highest  layers 
at  most  latitudes.  It  appears  that  the  use  of  the  alternative  data  in  the 
01  analyses  resulted  in  better  overall  analyses  at  the  lowest  and  highest 
layers  and  somewhat  poorer  analyses  in  the  mid-layers.  Since  the 
differences  are  most  evident  in  conventional  data  sparse  regions,  the 
additional  moisture  information  appears  to  be  an  improvement  to  the  first 
guess  at  the  bottom  and  top  layers,  and  leads  to  a  worsening  of  the  first 
guess  in  the  mid-layers.  These  conclusions  must  be  tempered  somewhat  by  the 
availabilty  of  rawinsonde  comparison  data  as  shown  in  Fig.  16,  and  by  the 
limited  sample  (10  cycles)  of  the  assimilation  period. 
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Pig.  20(a,b)  Zonal,  Tina  Average  Relativa  Humidity  for  (a)  ASAP  (Opt.  2) 
Analyses  and  (b)  Difference  between  (a)  and  Basis 
Observations  (Pig.  17(a)).  Contour  Interval  for  (a)  is  5%, 
for  (b)  is  2V  Dashed  Contours  Indicate  legative  Values. 
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Ill  an  attempt  to  assess  the  effect  of  the  alternative  relative  humidity 
information  on  an  individual  analysis.  Pigs.  21-23  were  prepared.  Pig.  21 
shows  the  locations  for  the  Vorthern  Hemisphere  of  the  observations  (a) 
which  survived  the  gross  and  buddy  checks,  and  (b)  which  are  actually  used 
in  the  analysis  for  2/17/79  12Z  for  (~860  ®b> .  Notice  once  again 
the  distinct  preference  for  the  use  of  rawinsonde  information  in  the 
analysis.  This  results  in  rawinsonde  rich  regions,  such  as  North  America, 
Europe,  and  northern  Asia,  which  are  left  exclusively  with  rawinsonde 
relative  humidities  with  which  to  perform  the  analysis.  Therefore,  we  would 
expect  the  effect  of.  the  corrections  on  the  first  guess  from  moisture 
options  1  and  2  to  look  quite  similar  over  these  regions.  This  is  confirmed 
by  comparing  Pig.  22(a)  (analysis  corrections  for  option  1)  with  Pig.  22(b) 
(analysis  corrections  for  option  2).  Note  the  relatively  similar  analysis 
corrections  over  the  three  aforementioned  rawinsonde  rich  areas.  Pig. 
23(a,b),  which  shows  the  corresponding  analysis  errors  for  the  two  moisture 
options,  indicates  by  their  similarly  low  values  of  around  8  percent  in 
these  regions  that  the  low  OESO  values  for  rawinsondes  have  led  to  low 
analysis  errors  throughout  the  assimilation  period  for  both  moisture 
options.  As  expected,  the  .primary  differences  between  the  two  sets  of 
analysis  corrections  should  occur  in  the  rawinsonde  sparse  areas.  Por 
example,  consider  the  region  bounded  by  latitudes  30*N  and  60*N  and 
longitudes  180*B  and  240*B.  This  region  of  the  northern  Pacific  had 
primarily  rawinsondes  used  in  the  analysis  on  its  northern  and  eastern 
edges,  but  predominantly  SWO  and  3DN  pseudo-observations  in  the  rest  of  the 
area  (see  Pig.  21(b)].  Notice  that  in  Pig.  22(b)  (moisture  option  2)  the 
analysis  corrections  are  smaller  for  the  predominantly  SWO,  3DN  part  of  the 
region  than  for  the  nearby  rawinsonde  rich  regions.  Because  of  the  complete 
lack  of  information  in  this  area  for  moisture  option  1,  the  analysis 
corrections  are  essentially  zero  in  Pig.  22(a).  As  expected  by  the 
difference  in  assigned  OBSD's  for  the  rawinsonde  and  alternative  relative 
tumidity  values,  the  corrections  to  the  first  guess  in  regions  dominated  by 
the  rawinsonde  are  significantly  different  from  the  corrections  in 
alternative  data  rich  regions,  with  the  corrections  in  the  rawinsonde 
dosdnated  areas  being  generally  larger.  A  better  measure  of  the  relative 
effect  is  depicted  in  Pig.  23  (a,b).  Here  it  is  seen  that  whereas  in 
rawinsonde  rich  areas  the  analysis  error  is  similar  for  the  two  analyses,  in 
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fi*.  21(a,b)  Location  of  (a)  Chackad  and  <b)  Oaad  Balativa  Humidity 

Obaarvations  for  ASAP  (Opt.  2)  Analysis  of  2/17/79,  12Z  at 
02  (-860  mb) .  R  •  Rawinsonda,  X  »  Surfaca  Uaathar 
Obsarvation  Infarrad,  •  *  3DRKPH  Inf air ad. 
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Pig.  22(a,b)  Analysis  Corrsctions  for  ASAP  Rslativs  Humidity  Analyses  (a) 
Opt.  1  and  (b)  Opt.  2  for  2/17/79,  12Z  at  a2  (~860 
mb) .  Contour  Interval  is  5%.  Dashed  Contours  Indicate 
Vegative  Corrections. 
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rawinsonde  sparse  but  alternative  data  rich  regions  the  analysis  error  is 
noticeably  reduced  below  the  maximum  allowed  (48  percent),  which  occurs  in 
all  of  the  rawinsonde  void  regions  remaining  essentially  uncorrected  through 
the  assimilation.  Thus,  it  appears  that  the  use  of  the  alternative 
information  leads  to  reduced  analysis  errors  over  no  analysis,  but  that  the 
corrections  in  areas  dominated  by  alternative  information  are  generally 
smaller  than  in  rawinsonde  rich  areas.  The  smaller  corrections  are  to  be 
expected  because  the  analysis  error  in  alternative  rich  regions  (typically 
20-24  percent)  is  larger  than  in  rawinsonde  rich  regions  (8  percent);  thus 
any  alternative  observations  would  have  corresponding  less  weight  in  the 
analysis.  Larger  OESD's  lead  to  less  weight  for  the  observations  (and  thus 
generally  a  smaller  total  correction) ,  which  in  turn  leads  to  higher 
analysis  errors.  Votice  in  Fig.  23(b)  that  in  all  of  the  alternative 
humidity  rich  areas,  the  analysis  error  (and  thus  the  FESD)  is  maintained  at 
a  level  approximating  the  OESD  for  850  mb  SWO's  and  3DH*s  (see  Table  13). 
Comparing  the  original  values  of  relative  humidity  FESD  used  at  the 
beginning  of  the  assimilation  (see  Table  4) ,  it  appears  that  the  alternative 
sources  of  humidity  information  are  effective  in  holding  the  FESD  at  a 
relatively  level  value  through. the  assimilation  for  areas  where  they  are 
relatively  plentiful. 

Because  the  3DVEPH  cloud  amounts  represent  an  analysis  themselves,  it 
is  not  necessary  to  perform  an  additional  analysis  to  move  them  to  the  grid 
locations.  The  horizontally  averaged  cloud  amount  are  available  on  a  50  km 
(quarter-mesh)  grid,  and  as  such  very  little  error  is  introduced  in  simply 
interpolating  bilinearly  the  gridded  cloud  amounts  (converted  to  cos 
((w/100)  •  X  cloud  amount])  to  the  analysis  points  at  each  a  layer. 

This  was  done  for  the  cloudy  regions  in  moisture  option  3,  and  the  relative 
humidity  inferred  from  these  cloudy  analysis  point  cloud  amounts  was  used  to 
directly  replace  the  first  guess  relative  humidity  at  each  analysis  time. 

At  non-cloudy  analysis  points,  an  OX  analysis  using  rawinsonde  relative 
humidity  only  was  performed  if  such  data  were  available  nearby,  otherwise 
the  first  guess  value  was  maintained.  Table  16  gives  the  percentage  of 
analysis  points  through  the  five  day  period  which  fell  into  each  of  the 
three  categories.  Since  the  3DH  replacement  occurred  at  cloudy  analysis 
points  only,  the  3DII  percentages  are  an  estimate  of  global  cloud  cover  at 


Table  16.  Percent  of  Analysis  Points  in  ASAP3 


(3DM)  Replacement  by  3 DR  Values  Occurred,  (01)  Rawinsonde 
01  Occurred,  (FG)  First  Guess  was  Maintained 


MB 

3DN 

275 

22 

337 

22 

436 

18 

574 

31 

724 

51 

862 

51 

962 

48 

01 

FG 

49 

29 

50 

28 

54 

28 

46 

23 

33 

16 

36 

13 

38 

14 

93 


each  0  layer.  Fig.  24(a-c)  shows  these  same  percentages  as  zonal  averages 
for  the  five  day  period.  According  to  Fig.  24(a),  cloud  cover  is  generally 
a  maximum  at  the  lower  levels  through  most  of  the  summer  hemisphere,  and  at 
the  724  mb  layer  in  the  higher  latitudes  of  the  winter  hemisphere.  At  the 
436  mb  layer  and  above,  cloud  cover  is  for  the  most  part  less  than  30 
percent,  and  it  is  at  these  layers  and  in  the  descending  Hadley  branch  that 
conventional  01  reaches  its  maximum  values  (Fig.  24b).  Maintenance  of  the 
first  guess  exceeds  50  percent  only  in  the  highest  three  layers  and  high 
latitudes  of  the  Southern  Hemisphere. 

Fig.  25(a,b)  shows  the  zonal,  time  average  of  the  five  day  option  3 
analysis  and  its  difference  from  observations  respectively.  Of  immediate 
interest  is  the  comparison  of  the  analysis-observation  differences  for 
moisture  option  2  (Fig.  20(b))  and  3  (Fig.  25(b)).  Notice  that  the  two 
moist  anomalies  in  the  ASAP2-OBS  figure  are  reduced  dramatically  in  the 
ASAP3-OBS  figure.  The  Southern  Hemisphere  anomaly,  centered  at  the  712  mb 
layer  at  30*S-40*S,  is  in  an  area  of  about  50  percent  replacement  by  the  3DM 
analysis  point  values  and  about  30  percent  conventionally  analyzed 
observation  sites  in  the  zonal  mean.  Thus,  it  appears  that  the  replacement 
3DN  values  had  an  improving  effect  on  the  humidity  analysis  at  this  locale. 
The  other  ASAP2-OBS  anomaly,  in  the  descending  Hadley  branch,  is  in  a  region 
of  about  20  percent  3DN  replacement  and  60  percent  conventional  01.  The 
anomaly  is  reduced  by  about  half  in  magnitude,  and  this  must  be  due 
primarily  to  the  conventional  01  without  the  use  of  the  alternative  data. 
This  would  reinforce  the  earlier  idea  that  this  anomaly  is  a  result  of  the 
worsening  of  the  analysis  by  the  use  of  the  alternative  information  used  in 
the  01  mode.  It  is  important  in  the  case  of  the  3DN  to  make  the  distinction 
between  use  as  pseudo-observations  in  the  01  as  opposed  to  its  use  in  a 
direct  replacement  of  the  first  guess  value.  In  the  former  case,  it  takes 
the  form  of  a  correction  whose  weight  is  roughly  equal  to  the  first  guess 
field  it  seeks  to  correct.  In  the  replacement  case,  it  completely  replaces 
the  first  guess,  acting  as  having  a  weight  of  1.0  vs.  the  first  guess  weight 
of  0.0.  Thus,  in  cloudy  regions  the  analysis  field  is  completely  dominated 
by  the  features  of  relative  humidity  inferred  from  cloud  amounts.  In  such 
areas  (indicated  by  over  50  percent  3DN  replacement  in  Fig.  25(a),  such  as 
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Fig.  24(a-c)  Zonal  Averages  of  Percent  Frequency  of  (a)  3DHEPH  Replacement, 
(b)  Conventional  01  Using  Rawinsonde  Observations,  and  (c) 
Maintaining  the  First  Guess  during  the  Ten  00.12Z  Cycles  of  the 
ASAP  (Opt.  3)  Relative  Humidity  Assimilation  Experiment. 
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Pig.  25 (a, b)  Zonal,  Time  Average  Relative  Humidity  for  (a)  ASAP  (Opt.  3) 

Analyses  and  (b)  Difference  between  (a)  and  Basis  Observations 
(Fig.  17(a)).  Contour  Interval  for  (a)  is  5%,  for  (b)  is  2%. 
Dashed  Contours  Indicate  Negative  Values. 


at  and  below  throughout  the  Southern  Hemisphere  and  at  a2~a4 
in  the  Northern  Hemisphere  extratropics),  the  3DH  replacement  performs  quite 
favorably  in  comparison  to  the  option  1  analyses  (Fig.  17(c)).  In  fact,  the 
moist  bias  occurring  at  the  lowest  layer  in  the  ASAP1-OBS  figure  is  reduced 
significantly  through  the  Southern  Hemisphere.  Yet,  the  mid-layers  in  the 
Northern  Hemisphere  extratropics  are  only  slightly  drier  (~2  percent)  in 
the  ASAP3  analyses  than  in  the  ASAP1  analyses.  Thus,  it  appears  that  the 
3DN  inferred  humidities  are  contributing  positively  to  the  analysis  when 
used  in  the  replacement  mode,  even  though  they  may  not  be  doing  so  when  used 
in  conjunction  with  SWO  inferred  humidities  in  the  01  mode.  In  the  latter 
case,  it  could  be  that  their  relatively  high  assigned  OESD  values  result  in 
a  damping  of  their  contribution,  being  masked  by  the  first  guess  which  is 
considered  equally  as  valid  by  the  analysis  procedure.  Looking  back  at 
Tables  13  and  14,  it  is  clear  that  if  the  OESD  for  rawinsonde  relative, 
humidities  were  raised  to  IS  percent  at  all  levels,  the  OESD  values  for  the 
alternative  observations  would  drop  3-5  percent,  which  would  give  them  more 
impact  in  the  analysis.  However,  the  fact  that  their  straight  replacement 
in  the  analysis  appears  to  be  making  a  positive  contribution  to  the 
assimilation  analyses  would  suggest  that  they  should  be  used  with  their 
fullest  effect,  namely  as  replacements  rather  than  observations  for 
interpolation.  Then  it  would  remain  to  be  seen  if  SWO  observations  by 
themselves  result  in  an  improvement  of  the  first  guess  in  the  01  mode. 

In  this  vein,  three  further  experiments  were  conducted  to  examine  the 
effects  of  3DN  replacement  and  SWO  01  analysis  in  areas  where  the  resulting 
analyses  can  be  compared  to  a  plentiful  supply  of  high  quality  observations, 
namely  rawinsondes.  It  would  then  be  assumed  that  alternative  information 
would  have  the  same  effect  on  the  analyses  in  conventionally  data  sparse 
areas.  Moisture  option  4  involved  the  use  of  3DN  replacement  in  cloudy 
regions  only  (i.e,  no  conventional  01  in  non-cloudy  regions).  Moisture 
option  5  allowed  no  correction  to  the  first  guess  anywhere  during  the  course 
of  the  assimilation  (to  act  as  an  experimental  control).  Moisture  option  6 
allowed  the  use  of  only  surface  weather  observations  in  the  01,  with 
maintenance  of  the  first  guess  elsewhere.  In  the  option  6  case,  the  OESD's 
for  SWO' s  were  reduced  to  a  level  corresponding  to  OESD's  of  15  percent  RH 
for  rawinsondes.  Thus  their  values  at  levels  1000-300  mb  were  assigned  as 


8.0,  16.0,  20.0,  19.0,  18.0,  16.0.  This  was  done  to  allow  them  to  have 
their  fullest  effect  in  the  analysis  in  order  to  better  assess  their  impact 
on  the  first  guess.  ^ 

RMS  differences  and  bias  against  the  basis  set  of  rawinsonde  relative  . 
humidity  at  00Z,  12Z  times  were  computed  and  plotted  in  Figs.  26  and  27. 

The  solid  curve  in  each  plot  (no  humidity  analysis)  represents  the  tendency 
of  the  forecast  model  in  its  development  of  relative  humidity  vs.  observed 
humidity.  In  the  RMS  curves,  there  is  a  clear  growth  of  RMS  difference  for 
the  "no  analysis"  curves,  especially  at  a  layers  3  and  above.  The 
corresponding  bias  curves  show  a  tendency  of  the  model  to  generate  slightly 
too  moist  a  forecast  at  o^,  too  dry  at  o^,  <>3,  and  substantially 
too  moist  at  ®.-o, . 

As  for  the  ASAP4  (3DN  replacement)  and  ASAP6  (SWO  01  only),  it  is  fair 
only  to  compare  each  against  the  "no  analysis"  curves,  and  not  against  each 
other.  This  is  because  the  ASAP4  analyses  represent  3DV  replacment  only 
where  clouds  occur,  which  would  show  no  geographic  preference  for  where 
comparative  rawinsondes  are  located.  Since  the  ASAP4  curves  represent  the 
effect  of  changes  from  the  background  field  only  in  the  area  of  comparative 
rawinsondes,  this  would  represent  in  general  only  about  half  of  the  total 
effect  of  3DN  replacement  globally.  On  the  other  hand,  a  majority  of  the 
SWO's  occur  in  areas  rich  in  rawinsondes,  especially  over  certain  land 
regions  (see  for  example,  Fig.  19(a)).  Therefore,  most  of  their  effect  on 
the  analysis  will  be  felt  in  areas  Where  there  are  rawinsondes  against  which 
to  measure  their  impact  on  the  analysis. 

Looking  first  at  the  RMS  and  bias  for  the  3DM  replacement  case,  the  RMS 
curves  appear  to  represent  an  improvement  to  the  uncorrected  first  guess 
("no  analysis"  curves)  at  all  layers  except  possibly  where  on  the 
whole  it  shows  no  clear  improvement  over  the  model  generated  humidity 
field.  The  bias  curves  for  3DN  replacement  show  a  lessening  of  the  bias 
from  the  "no  analysis"  case  at  and  a^-Cy,  and  an  introduction  of  a 
substantial  moist  bias  at  a 2,  a^.  At  these  two  layers,  the  tendency 
of  a  model  dry  bias  is  overcompensated  for  by  an  apparently  larger  moist  bias 
due  to  the  3 DR  replacement.  That  the  3DR  replacement  curves  show  a  moist 
bias  relative  to  the  model  generated  humidity  fields  is  somewhat  surprising 
in  light  of  Table  12,  which  shows  a  dry  bias  at  all  levels  in  the 
3DR-rawinsonde  collocation  study  using  the  Tibaldi  method.  However,  the 
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Fig.  26  SMS  Differences  on  a  Layers  between  Relative  Humidity  Fields 
and  Basis  Observations 
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more  important  matter  here  is  how  the  forecast  model  responds  to  the 
replacements.  That  is,  does  the  3DM  replacement  lead  to  a  more  moist  6h 
forecast  at  these  layers  than  would  occur  if  no  replacement  took  place?  To 
determine  the  relative  bias  of  the  3DN  replacement  assimilation  first  guess 
fields,  they  were  compared  with  the  same  basis  rawinsonde  relative 
humidities  at  all  00,  12Z  times,  and  their  bias  is  plotted  on  the  c 2, 
o,  bias  plots.  The  results  make  it  clear  that  at  these  two  layers,  the 
3DM  replacement  analyses  lead  to  6h  forecast  fields  that  are  more  moist  than 
an  uncorrected  first  guess.  In  the  GSM,  the  modest  dry  bias  in  the 
uncorrected  first  guess  is  replaced  by  a  more  substantial  moist  bias  in  the 
6h  forecast  when  the  3DN  humidities  from  cloudy  areas  are  used  in  place  of 
the  model  generated  first  guess  for  the  initial  conditions.  This  relative 
moist  bias  is  retained  in  later  analyses  as  the  replacement  method  is 
ineffective  in  removing  it.  At  all  other  layers  where  the  uncorrected  first 
guess  was  too  moist,  the  3DN  replacement  analyses  were  definitely  and 
consistently  less  moist  than  their  counterpart  first  guess,  so  in  these 
layers  the  replacement  was  effective  in  reducing  the  model  moist  bias. 

The  SWO  analyses,  in  considering  both  RMS  and  bias  at  each  layer,  are 
actually  in  poorer  agreement  with  rawinsonde  relative  humidity  than  the  "no 
analysis"  fields  at  all  layers.  Whereas  the  RMS  curves  for  SWO  show  better 
agreement  at  all  layers  except  a the  bias  curves  show  a 
greater  bias  at  all  layers  except  It  appears  that  the  use  of  SWO 

values  in  conventionally  data  sparse  regions  would  lead  to  an  artificially 
too  moist  depiction  of  the  atmospheric  state  except  for  possibly  at  a 2> 
where  the  RMS  fit  is  much  poorer  than  the  uncorrected  first  guess.  On  the 
basis  of  these  results,  one  would  be  better  off  to  leave  the  GSM  uncorrected 
in  such  areas  rather  than  use  the  SWO's  to  attempt  a  correction.  This 
conclusion  of  course  is  model  dependent,  and  cannot  necessarily  be  applied 
to  other  global  forecast  models  used  for  data  assimilation. 

Finally,  it  is  of  interest  to  examine  the  model  generated  humidity 
fields  in  an  averaged  sense.  Fig.  28(a)  is  a  zonal,  time  average  of  the 
uncorrected  (ASAP5)  6h  forecast  fields  averaged  over  all  00,  12Z  times  in 
the  five  day  assimilation.  Fig.  28(b)  represents  the  difference  between 
this  plot  and  the  same  plot  for  the  basis  rawinsondes.  Note  the  tendency 
for  the  model  forecast  to  be  too  moist  in  the  lowest  and  highest  layers.  A 
comparison  of  Fig.  28(b)  with  Fig.  17(b)  (ASAP1)  and  Fig.  25(b)  (ASAP3) 


reveals  that  a  good  analysis  scheme  will  on  the  average  lead  to  a  good  time 
averaged  depiction  of  the  relative  humidity  when  produced  by  an  objective 
data  assimilation  system. 

How  that  it  is  acknowledged  that  certain  types  of  humidity  information 
can  lead  to  good  analyses  of  relative  humidity  on  the  average,  the  question 
of  whether  this  in  turn  leads  to  a  better  long  term  forecast  of  relative 
humidity  must  be  addressed.  At  this  point,  the  particular  characteristics 
of  the  model  being  used  to  produce  such  a  forecast  become  the  most  important 
factor  in  determining  the  answer  to  the  question.  Can  a  somewhat  coarse 
resolution  model  with  modest  physical  parameterizations  such  as  the  GSM  be 
sufficiently  sensitive  to  an  initial  relative  humidity  field  so  as  to 
produce  a  more  highly  verified  forecast  from  an  initial  field  that  more 
closely  agrees  with  its  objective  verification  data?  Or.  are  the  initial 
mass  and  motion  fields  more  important  in  determining  the  forecast 
distribution  of  humidity?  In  an  attempt  to  answer  these  questions,  the  48h 
humidity  forecast  results  at  12h  intervals  based  on  the  2/17/79  12Z  analysis 
were  compared  with  basis  rawinsondes  and  RMS  and  bias  curves  were  produced 
and  plotted  in  Figs.  29  and  30.  Fig.  31  is  a  plot  of  the  RMS  differences 
and  bias  for  the  initial  RH  analyses  at  2/17/79  12Z  for  ASAP  and  F3A  along 
with  the  RMS  and  bias  for  all  ten  00,  12Z  times.  It  is  clear  from  both  the 
RMS  curves  and  bias  curves  that  the  initially  closer  fit  of  the  ASAP1 
analysis  to  the  observations  than  that  for  F3A  as  seen  in  Fig.  31  did  not  in 
this  case  lead  to  a  closer  fit  for  the  forecast  in  Figs.  29  and  30. 

Hotice  first  in  the  RMS  curves,  by  comparing  the  12h  forecast  curves 
with  the  initial  curves,  that  the  forecast  model  quickly  moves  the  RMS 
values  to  the  20-30  percent  range  where  the  initial  F3A  values  were.  In 
fact,  after  12h  the  RMS  differences  for  F3A  based  forecast  are  slightly 
lower  at  all  layers  except  the  lowest,  and  continues  to  be  lower  even  after 
48h.  Clearly,  the  better  initial  RMS  fit  in  the  ASAP1  analysis  had  no 
positive  impact  upon  the  resulting  forecast,  and  the  fact  that  the  ASAP1 
forecast  RMS  values  are  actually  greater  than  those  for  the  F3A-based 
forecasts  suggests  that  other  differences  (mass-motion)  between  the  two 

analyses  had  a  greater  effect  than  initial  humidity  on  the  humidity 

% 

forecast.  The  three  ASAP  based  forecasts  show  little  difference  in  RMS; 
however,  the  ASAP2  does  appear  to  have  higher  values  than  the  other  two. 


Differences  on  a  Layers  between  Forecast  Relative  Humidity  Fields  and  Basis  Observations 


Observation  Bias  for  o  Layer  Relative  Humidity  (Basis  Observations) 
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SMS  Differences  between  Analysis  Fields  and  Observations  and  Field-Observation  Bias  for 
Layer  Relative  Hunidity 


A  look  at  the  forecast  bias  for  the  F3A  based  forecast  as  compared  to 
the  ASAP  based  forecasts  in  Fig.  30  reveals  that  the  F3A  based  forecasts 
tend  to  be  less  dry  in  the  four  lowest  layers  and  less  moist  in  the  three 
highest  layers  after  12h.  By  48h,  the  F3A  and  ASAP  based  forecasts  all  show 
about  the  same  amount  of  dry  bias  in  the  lowest  four  layers,  but  the  F3A 
based  forecast  shows  less  of  a  moist  bias  in  the  upper  three  layers.  The 
reason  for  these  trends  is  clear  from  a  comparison  with  the  initial  biases 
in  Fig.  31.  The  F3A  analysis  has  a  clear  moist  bias  in  the  lowest  three 
layers  where  the  forecast  model  tends  to  generate  humidities  that  are  too 
low,  while  in  the  next  three  layers  the  F3A  analysis  shows  a  pronounced  dry 
bias  where  the  model  tends  to  generate  humidities  that  are  too  high.  Thus, 
when  the  model  moves  the  bias  to  the  left  (drier)  in  the  lowest  three  layers 
and  to  the  right  (more  moist)  in  the  next  three  layers,  the  F3A. based 
forecast  shows  lower  bias  because  its  initial  fields  compensated  in  sign 
more  greatly  than  did  the  lower  bias  ASAP1  fields.  It  is  evident  from  these 
results  that  the  model  establishes  its  own  bias  profile  configuration 
without  regard  to  the  configuration  of  the  initial  bias.  Here  again,  very 
little  difference  is  seen  among  the  biases  for  the  three  ASAP  based 
forecasts,  although  the  ASAP2,  ASAP3  biases  are  higher  than  ASAP1  after  48h. 

It  is  apparent  in  this  case  that  for  the  GSM  under  these  circumstances, 
closeness  of  fit  of  the  initial  humidities  with  observations  did  not  result 


in  better  verification  with  observations  in  the  ensuing  forecast. 

Presumably  the  GSM  in  its  present  configuration  is  not  capable  of  preserving 
the  better  initial  verification  of  the  ASAP  humidity  fields  in  its 
forecast.  It  may  be  that  the  difference  between  the  F3A  and  ASAP  mass  and 
motion  fields  plays  a  bigger  role  in  determining  the  forecast  distribution 
of  humidity  than  does  the  initial  humidity  distribution.  This  was  a 
conclusion  reached  by  Lejenas  (1979)  in  a  study  using  a  regional  model  for 
precipitation  forecasting.  He  concluded  that  vertical  velocities  in  the 


initial  wind  field  are  mo^e  important  than  fine  detail  in  the  moisture  field 
for  short  term  (<24h)  forecasts,  and  that  for  global  long  term  forecasts 
initial  humidity  specification  is  of  little  importance  to  the  forecast 
humidity.  Similar  conclusions  were  reached  by  Smagorinsky  et  al.  (1970), 
who  found  that  after  about  6  hours,  there  were  generally  no  significant 
differences  in  forecast  precipitation  rate  between  forecasts  involving 
various  moisture  specifications.  However,  in  the  present  study  cloud 


mmwm 


distribution  is  of  more  importance,  which  in  turn  suggests  that  humidity 
rather  than  precipitation  is  the  more  important  factor.  In  response  to  the 
conclusion  by  Lejenas,  it  must  be  stated  that  the  degree  to  which  that 
conclusion  is  true  is  dependent  upon  the  forecast  model,  as  was  stated 
previously.  To  determine  whether  the  mass-motion  fields  might  have  the 
greater  effect  on  the  humidity  forecast  in  the  case  of  the  GSM,  the  relative 
humidity  fields  in  the  F3A  and  ASAP1  analyses  were  switched  and  new  48h 
forecasts  were  generated.  Forecast  SMS  difference  and  bias  in  comparison  to 
basis  observations  were  generated  and  are  displayed  in  Figs.  32  and  33.  To 
distinguish  more  clearly  the  effect  of  switching  the  relative  humidity 
between  the  two  initial  analyses  on  the  resulting  forecasts,  the  RMS  and 
bias  curves  from  Figs.  29  and  30  are  included. 

Considering  first  the  RMS  curves,  it  is  apparent  in  the  comparison  of 
the  two  dashed  (F3A  mass-motion)  curves  that  using  the  better  fitting 
initial  ASAP  humidities  had  very  little  effect  on  the  humidity  forecast. 

The  relative  humility  forecast  from  the  F3A  analysis  improved  only  slightly 
(at  most  levels  less  than  1  percent  RMS)  when  the  better  fitting  initial 
ASAP1  humidity  analysis  replaced  the  F3A  humidity  analysis.  However,  the 
relative  humidity  forecast  from  the  ASAP  analysis  became  more  substantially 
worse  (at  most  layers  more  than  2  percent  RMS)  when  the  poorer  fitting 
initial  F3A  humidity  analysis  replaced  the  initial  ASAP1  analysis.  This 
suggests  that  closeness  of  fit  of  the  original  humidity  analysis  played  a 
greater  role  in  the  humidity  forecast  for  the  ASAP  based  forecast  than  for 
the  F3A  based  forecast.  In  general,  the  humidity  forecasts  from  the  F3A 
height-wind  fields  have  a  lower  RMS  than  the  forecasts  from  the  ASAP1  fields 
regardless  of  which  initial  relative  humidity  analysis  was  used.  This 
indicates  that  the  initial  mass-motion  analysis  has  a  greater  role  than 
initial  relative  humidity  in  determining  GSM  relative  humidity  forecast 
accuracy.  However,  the  fact  that  the  use  of  ASAP1  initial  humidities  in  the 
ASAP1  analyses  led  to  better  humidity  forecasts  than  using  F3A  initial 
humidities  indicates  that  accuracy  in  initial  humidity  analyses  can  play  a 
secondary  role  in  the  resulting  humidity  forecast. 

The  bias  curves  in  Fig.  33  confirm  the  earlier  conclusion  that  the  GSM 
establishes  its  own  moisture  bias  regardless  of  initial  humidity  bias 
configuration.  The  comparisons  between  the  two  F3A  based  forecasts  (dashed 
curves)  and  two  ASAP1  based  forecasts  (solid  curves)  are  shown  only  for  36h 
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Difference*  on  0  Layer*  between  Forecast  Relative  Humidity  Field*  and  Basis  Observation* 


Forecast  Field  -  Observation  Bias  for  o  Layer  Relative  Husiidity  (Basis  Observations) 


and  48h  when  the  model  has  clearly  established  its  bias  configuration,  in 
both  pairs  of  curves,  the  lower  dry  biases  and  upper  moist  biases  are  lower 
for  the  forecasts  which  use  the  initial  F3 A  relative  humidities.  This  is 
due  to  the  earlier  stated  fact  that  the  initial  F3A  bias  configuration  more 
substantially  compensated  for  the  model  generated  biases.  Here  again,  it  is 
clear  that  the  mechanics  of  the  model  overwhelmingly  determine  the  bias  of 
the  forecast  humidity  fields,  and  the  initial  humidity  biases  play  very 
little  role  in  affecting  the  trends  (as  opposed  to  actual  values)  of  the 
model  generated  bias.  To  determine  the  role  of  the  initial  mass-motion 
analysis  in  the  forecast  humidity  bias,  we  compare  the  two  curves  with  F3A 

relative  humidity  ( -  and  x - x)  on  one  hand,  and  the  two  curves  with 

ASAP1  relative  humidity  ( -  and  x - x)  on  the  other.  In  both  cases, 

the  solid  curves  (ASAF1  mass-motion)  show  moderately  less  bias  at  36h  than 
their  F3A  mass-motion  counterparts  at  all  layers  except  the  top.  However, 
much  of  the  difference  is  lost  by  48h  and  in  some  cases  is  actually 
reversed.  At  most,  the  mass-motion  analyses  play  a  secondary  role  in 
determining  GSM  humidity  forecast  bias,  with  the  model  itself  playing  the 
greatest  role. 

VII.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  STUDY 

On  the  basis  of  the  results  of  the  five  day  assimilation  experiments 
described  in  this  report,  the  following  conclusions  are  given. 

1)  Observation  error  standard  deviation  chosen  for  mass  and  wind  fields 
plays  a  crucial  role  in  determining  degree  of  fit  of  the  analysis  to 
observations.  A  good  index  of  relative  observation  error  influence  is  the 
ratio  of  normalized  observation  error  for  height  to  the  same  for  winds. 

2)  The  baseline  AFGL  GSM  used  in  the  data  assimilation  system  produces 
a  large  positive  height  bias  at  the  a  *  0.05  level,  which  is  a  result  of  a 
warm  temperature  bias  in  the  layer  topped  by  this  level.  Because  of  the  lack 
of  high  quality  height  and  wind  data  available  at  or  near  that  level,  the 
ASAP  analysis  is  unable  to  correct  this  bias  substantially  at  each  analysis 
time.  As  a  result,  the  bias  grows  in  time  during  the  course  of  the 
assimilation.  However,  the  unabated  growth  of  the  bias  in  the  assimilation 
analyses  is  restricted  to  the  layer  above  o  =»  0.10. 

3)  Initial  fit  to  observations  appears  to  play  a  greater  role  in 
forecast  verification  for  winds  than  it  does  for  heights  in  48h  forecasts 


produced  by  the  GSM.  A  close  fit  to  heights  in  the  initial  conditions  does 
not  necessarily  result  in  a  better  height  forecast  verification,  whereas  the 
GSM  does  seem  to  preserve  the  relationship  between  the  relative 
observational  data  fits  for  winds  for  several  initial  conditions  through  the 
forecast  period. 

4)  Fit  of  univariate  relative  humidity  analyses  to  observations  is  very 
much  dependent  upon  the  observation  error  standard  deviations  assigned  to 
the  various  observation  types.  Small  observation  error  standard  deviation 
(OBSD)  values  for  a  certain  data  type  result  in  a  closer  fit  of  the  analysis 
to  that  type  of  observational  data. 

5)  The  relatively  low  values  of  OESD  for  rawinsonde  relative  humidity 
observations  and  the  relatively  high  values  of  OBSD  for  alternative  humidity 
sources  lead  to  a  strong  preference  for  use  of  rawinsonde  observations  in  a 
mixed  analysis.  Thus,  in  areas  with  even  a  few,  somewhat  scattered 
rawinsonde  observations,  the  analysis  procedure  will  use  those  before  using 
closer  and  more  densely  concentrated  alternative  humidity  estimates.  In 
regions  rich  in  rawinsonde  observations,  alternative  humidity  sources  are 
effectively  excluded  from  use  in  the  analysis. 

6)  Though  the  magnitude  of  relative  humidity  corrections  may  be 
appreciably  smaller  in  alternative  humidity  rich  areas  than  in  their 
rawinsonde  rich  counterparts,  the  analysis  error  in  the  alternative  data 
rich  area  is  appreciably  reduced  from  its  value  when  no  corrections  are  made 
in  the  area.  Thus,  the  analysis  error  in  alternative  data  rich  areas  is  at 
a  level  above  that  for  rawinsonde  rich  areas  but  significantly  below 
climatological  limits  of  the  assimilation  system. 

7)  3DVEPH  replacement  of  first  guess  relative  humidity  has  a  more 
positive  effect  overall  on  the  analysis  than  the  use  of  3DHEPH  and  surface 
weather  observation  (SWO)  inferred  humidities  used  as  pseudo-observations  in 
an  analysis  mode.  This  suggests  that  since  the  3D  nephanalysis  is  already  a 
form  of  humidity  analysis,  3DNEPH  inferred  humidities  should  be  used  in  a 
replacement  mode  rather  than  in  an  analysis  correction  mode  to  achieve  their 
maximum  positive  effect. 

8)  Use  of  3DNEPH  inferred  relative  humidities  in  cloudy  regions  as  a 
direct  replacement  of  first  guess  humidities  represents  an  improvement  over 
no  humidity  correction  to  the  first  guess  at  layers  and  a ^  - 

<ty  At  Og  and  °3 •  use  of  3DMEPH  humidities  as  first  guess 
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replacements  leads  to  a  significant  moist  bias  in  the  ensuing  6h  GSM 
forecast,  where  a  modest  dry  bias  had  occurred  when  no  replacement  had  taken 
place.  The  replacements  at  these  layers  were  not  effective  in  reducing  the 
moist  bias  in  later  analyses. 

9)  Use  of  surface  weather  observation  inferred  humidities  in  an  01 
analysis  mode  represents  a  worse  analysis  overall  than  no  humidity 
correction  to  the  GSM  first  guess.  This  is  due  to  the  introduction  of  a 
significant  moist  bias  over  the  "no  analysis”  case  at  all  layers  except 
<t2,  where  the  RMS  differences  with  rawinsonde  relative  humidity  were 
excessive.  This  conclusion  is  dependent  upon  the  algorithm  used  to  infer 
humidities  from  SWO's  and  on  the  particular  assimilation  model.  In  the  case 
of  this  study,  the  Tibaldi  (ECMWF)  algorithm  is  used  for  the  former  and  the 
AFGL  GSM  (physics  as  in  Sela,  1980)  is  used  for  the  latter. 

10)  A  closer  fit  to  observations  in  one  initial  relative  humidity 
analysis  over  another  analysis  does  not  necessarily  result  in  a  better 
verification  against  observations  of  a  48h  forecast  when  the  GSM  is 
initialized  by  the  two  humidity  analyses.  Initial  RMS  fit  to  data  in  the 
relative  humidity  analysis  has  less  of  an  effect  on  the  humidity  forecast 
RMS  fit  than  does  the  initial  mass-motion  field.  The  GSM  tends  to  create 
its  own  humidity  forecast  bias  configuration  (too  dry  in  the  lower  layers, 
too  moist  in  the  upper  layers)  without  the  initial  humidity  bias  having  any 
effect  on  the  production  of  these  bias  trends.  Unfortunately,  it  was 
discovered  after  the  completion  of  the  experiments  that  the  version  of  the 
GSM  used  both  in  assimilation  and  forecasting  had  effectively  no  moist 
convective  adjustment.  A  companion  journal  article  (submitted  to  Monthly 
Weather  Review)  includes  additional  results  for  cases  when  the  moist 
convective  adjustment  scheme  is  utilized  as  described  by  Sela  (1980).  In 
general,  the  present  version  of  the  GSM  is  generally  insensitive  to 
specification  of  initial  relative  humidity  conditions  in  its  generation  of 
humidity  forecasts.  It  may  be  that  later  upgrades  in  the  GSM  moist 
convection  parameterization  (NMC  1983,  NMC  1985),  or  new  parameterizations 
of  moist  convection,  boundary  layer  physics,  and  radiation  developed  for 
ARGL  will  lead  to  an  improvement  in  model  sensitivity  to  initial  humidity 
distribution.  This  would  make  efforts  to  analyze  more  accurately  both 
conventional  data  rich  and  conventional  data  sparse  (alternative  data  rich) 
regions  worthwhile  in  that  they  would  lead  in  turn  to  more  accurate  humidity 


forecasts  in  both  types  of  regions. 

The  research  activities  described  in  this  report  are  a  first  attempt  to 
use  the  AFGL  GSM  as  a  global  data  assimilation  model.  As  such,  the  results 
of  this  study  answer  some  questions  about  data  assimilation  possibilities  at 
AFGL,  but  suggest  many  more  questions  that  need  to  be  answered  about  such  a 
system.  The  following  list  of  suggestions  for  further  study  is  by  no  means 
exhaustive,  but  merely  represents  the  most  obvious  starting  point  for 
further  exploration  of  the  role  of  various  types  of  observations  in 
four-dimensional  data  assimilation. 

1)  In  light  of  the  use  of  height  on  a  layers  as  the  mass  field 
analysis  medium,  more  study  of  the  best  means  of  converting  height  analysis 
corrections  to  tesq»erature  corrections  should  be  undertaken. 

2)  The  cause  of  the  large  warm  temperature  bias  at  a  layer 

A  A 

(i.e.,  the  layer  between  o  *  0.1  and  a  =  0.05)  in  the  12  layer  version 
of  the  GSM  should  be  explored  and  resolved. 

3)  At  least  two  more  days  of  assimilations  following  the  five  day 
experimental  period  described  in  this  study  should  be  carried  out  to  verify 
that  the  RMS  differences  between  analysis  winds  and  observational  winds  do 
not  grow  in  magnitude  beyond  their  fifth  day  level. 

4)  The  effect  of  global  average  temperature  changes  in  the  assimilation 
on  global  averaged  surface  pressure,  diagnosed  from  analyzed  a  layer 
temperatures ,  should  be  studied  in  more  detail  to  verify  the  realism  of  the 
derived  surface  pressures. 

5)  More  long  term  forecasts  should  be  run,  especially  from  initial 
fields  later  in  the  assimilation  period,  to  verify  conclusions  about 
forecast  characteristics  of  the  GSM. 

6)  A  new  collocation  study  of  FGCE  II-B  satellite  derived  layer 
precipitable  water  values  with  rawinsondes  should  be  conducted  using  the 
reprocessed  data  set.  A  In  q  vs.  In  p  assumption  should  be  used  in  each 
report  layer  to  convert  layer  precipitable  water  to  q. 

7)  Error  statistics  for  satellite  soundings  from  collocation  studies 
with  rawinsondes  for  both  Z  and  RH  should  be  conducted  for  land  and  ocean 
based  satellite  soundings  for  a  comparison  of  error  levels  with  those 
obtained  previously  for  all  soundings.  If  the  land  based  levels  are  found 
to  be  of  comparable  quality  to  ocean  based  soundings,  such  soundings  should 
be  included  in  the  assimilation  data  pool. 


8)  Other  methods  of  horizontal  and  vertical  compaction  of  3DEMPH  cloud 
amounts  should  be  tested  in  rawinsonde  humidity  collocation  studies  in  order 
to  maximize  the  amount  of  information  that  can  be  derived  from  the  cloud 
analysis  while  minimizing  the  collocation  differences  with  rawinsonde 
humidities.  Also,  the  collocation  study  should  be  carried  out  for  the  June 
period  to  verify  that  the  preference  for  the  Tibaldi  method  is  not 
seasonally  dependent. 

9)  More  long  term  humidity  forecasts  should  be  generated  using  the  GSM 
and  ASAP,  FGGE  III-A  initial  conditions,  to  better  ascertain  the  sensitivity 
of  the  model  to  initial  conditions  and  to  gain  a  better  understanding  of  its 
humidity  forecast  characteristics. 

10)  A  new  study  has  been  initiated  elsewhere  to  study  the  effect  of 
using  satellite  derived  data  (both  temperature  and  moisture)  in  the  ASAP  01 
procedures  on  the  data  assimilation  analyses  while  using  new  GSM 
parameterization  schemes  for  boundary  layer,  moist  convection,  and  radiation 
physics.  Once  these  studies  and  tests  have  been  completed  and  the  upper 
level  temperature  bias  corrected  in  the  GSM,  the  ASAP  experiments  ASAP1 
(conventional  RH  only),  ASAP5  (no  RH  analysis),  ASAP 4  (3DM  replacement 
only),  and  a  combination  of  conventional  RH  01  where  such  data  are  fairly 
plentiful  and  3 DM  replacement  where  possible  elsewhere,  should  be  run  in  an 
assimilation  mode.  Forecasts  generated  from  these  fields  should  be  examined 
to  see  if  humidity  forecasts  from  these  GSM  versions  benefit  from  better 
specification  of  initial  humidity  in  both  conventional  data  rich  and 
conventional  data  sparse  (alternative  data  rich)  regions. 


APPENDIX  A.  Flattery  Algorithm 


Given  the  system  of  linear  constraints  imposed  on  the  interface  temperatures 

A 

{Tj^} ,  k  *  1,  . ...  K  +  1  and  the  layer  temperatures  {T^}.  k  =  1,  ....  K.  We 
assume : 

(1)  that  the  layer  temperature  Tv  is  the  arithmetic  average  of  the  bounding 

*  A 

interface  temperatures  T^,  T^+1 

Tk  3  2  <Tk +  Tk+i}  •  k  ■ 1 . K;  (A1) 

(2)  that  temperature  in  the  layer  bounded  by  two  pressure  levels  p^, 

Pk+1  is  linear  in  In  p,  so  that  the  interface  temperature  T^+1  is  given 

in  terms  of  the  layer  temperatures  and  Tk+1  as 

T  =  T  W  +  T  W 

k+1  k  L,  k+1  k+1  U,  k+1 

where 


Hence 


Similarly 


Here,  the  subscript  k  designates  level  in  the  vertical  and  is  chosen  to 
increase  upward .  The  geometrical  configuration  may  be  depicted  as  follows: 


When  p^+1  ■  0,  we  cannot  employ  this  method.  There  are  different  alternatives 
here: 

*  I  *  m 

(1)  Introduce  pR+1  ,  which  is  between  p^  and  pR+1 

(2)  Assume  that  the  layer  above  pR  is  isothermal  so  that 


The  constraints  listed  above  may  be  put  into  matrix  form  given  by 


AT  a  BT 

* 

where  matrices  A,  B  and  vectors  T,  T  are  defined  as  follows: 


h  » H 


To  solve  Bq.  (A5)  in  accordance  with  the  least-squares  principle  we  proceed  as 
follows: 

(A)  When  {T^}  is  given, 

T  *  <ATA)-1ATBT  (A6) 

(B)  When  {X.  }  is  given, 

K  T  -1  T  * 

=  (B  B)  B  AT  . 


T 


(A7) 


APPENDIX  B.  Computation  of  Surface  Pressure  from  Analyzed 
Heights  of  a  Levels  and  Terrain  Height 


The  analysis  of  temperature  implies  a  correction  to  the  geopotential  height  at 

all  a  levels.  This  in  turn  implies  a  change  in  the  surface  pressure  at  the 

fixed  height  of  the  model  terrain.  To  estimate  the  update  surface  pressure  p 

*  * 

given  the  original  geopotential  ♦  and  update  geopotential  $  of  the  a 
levels,  consider  the  hydrostatic  equation  in  the  form 

3*/31np  =  -RT. 

Rewriting  in  the  form  pfr  =  -RT,  ^  noting  that 

ap 


a*  _  a±  a® 

3p  ~  3  a  3p 

where  a  =  p/p  and  3o/dp  *  1/p*,  we  have 


a*  a* 


3  a 


Thus,  —  =  -RT,  and  if  we  assume  T  varies  linearly  in  In  a,  then  #  must 

oin  a 

vary  quadratically  in  In  a: 

2 

♦(In  o)  =  C  +  A  In  o  +  B  In  a. 


If  we  denote  ♦. ,  ♦  ,  ♦  to  be  the  first  guess  geopotential  at  levels 

AAA  Z  J  A  * 

°2*  °3 ’  °4 ’  the  three  a  levels  immediately  above  =  1.0,  then  by  writing 
the  above  quadratic  expression  in  terms  of  these  three  levels,  and  solving 
the  resulting  equations  for  A,  B,  C,  we  have 


“*3  - 


♦2)/(ln  o3 


-  In 


°2>  -  <*4  - 


♦3)/(ln  o4 


-  In 


o3>J/(ln  o2 


-  In 


<*4> 
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0^  A  A  A  0k 

A  -  [<t4  -  *3)  -  B(ln  o4  -  In  o3)J/(ln  <»4  -  In  o3> 


C  ■  #3  -  A  In  o3  -  B  In  «3  . 


Similar  expressions  exist  for  A  ,  B  ,  C  in  terms  of  the  geopotentials 

•  *  2 

,  #  of  the  analyzed  field.  Using  these  values,  we  can  evaluate  the 

<3  ^ 

0k 

computation  terrain  geopotential  *A,  *A  at  *  1.0  from  both  fields: 


*  2*  c  ...  '  2  “ 

»  C  +  A  In  Oj  +  B  In  o3  »  C  ;  #A  =  C  +  A  In  +  B  In  «  C 

Given  the  first  guess  surface  pressure  pA  on  the  actual  model  terrain 
geopotential  #A,  we  can  estimate  the  computational  values  of  surface 
pressure  from 


pi  -  pA  exp[ (f A  -  t*)/(RT)l  ;  p®  -  pA  exp[(*A  -  )/(RT  )J 


where 


*  "  2  <T1  +  T<5)  *’  T°  "  <#2  ”  *£>/tR<ln  ai  ~  ln  <*2)] 


—  '  1  '  *  .  *  «  0k 

T  »  j  (Tx+  Tc  )  ;  Tc  -  (#2  -  )/[R(ln  -  ln  o2>J 


Then  the  final  estimate  of  updated  surface  pressure  is 


P*  *  P*  +  <p£  -  P*> 


t-p,*  ♦  I.' 


APPENDIX  C.  NMC  Subroutine  PTOSIG 


The  subroutine  generates  estimates  of  values  of  dynamic  variables  such  as 
temperature  T,  wind  speeds  u,  v,  and  specific  humidity  q,  all  of  which  are 
defined  as  layer  variables.  The  underlying  assumptions  are: 

(1)  The  geopotential  of  a  pressure  surface  is  a  quadratic  function  of  In  p. 

(2)  The  winds  and  natural  logarithm  (In)  of  specific  humidity  are  linear 
functions  of  In  p. 

(3)  The  temperature  T  in  the  layer  is  estimated  from  +  at  the  surrounding 
model  levels  using  the  hydrostatic  equation 

3^/3  In  p  *  -RT . 

A.  Calculation  of  Layer  Temperature 

The  first  step  is  to  calculate  the  geopotential  at  the  model  layer 
interfaces  (levels)  from  geopotential  42,  . . . ,  and 
temperature  ^ ,  T2 ,  . . . ,  T^  known  at  the  pressure  levels  p1 ,  p2 , 

....  p  .  Assuming  ♦  is  quadratic  in  In  p,  we  have 

la 


♦(In  p)  *  ♦  (In  p)  +  A(ln  p  -  In  p)  +  B  (In  p  -  In  p) 

2 


(Cl) 


for  p2  <  p  <  pj.  Here  A  and  B  are  constants  and 


In  p  -  ~  (In  p^  -f  In  p2) 


(C2) 


By  the  hydrostatic  assumption,  we  have 


d*/ d  In  p  -  -RT  *  A  ♦  B(ln  p  -  In  p) 
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a2+/aun  p)2  -  -R(aT/a  in  p)  -  b. 

The  last  relation  in  tha  second  aquation  states  that  T  oust  be  a  linear 
function  in  In  p;  hence 

dT/a  In  p  »  (Tj-Tjl/dn  Pj^  -  In  p2>  in  p2  <  p  <  pj 

and 

B  -  -  t2)/(ln  px  -  In  p2>.  (C3) 

■ext,  by  substituting  ♦^  +2  in  Bq.  (Cl)  we  obtain 

♦  *  ♦(In  p)  +  A (In  p  -  In  p)  +  B  (In  p  -  In  p)2  (C4) 

1  2  1 

♦  *  ♦(In  p)  +  A(ln  p  -  In  p)  +  B  (In  p  -  In  p)2.  (C5) 

2 

Adding  these  two  and  noting  the  definition  of  In  p  in  Eq.  (C2),  we  find 

1  (♦,  +  ♦,)  «  4>(ln  p)  +  B  {(In  p  -  1  (In  p  ♦  In  p  ]2  ♦ 

2  1  2  2  1  2 

[In  p.  -  1  (In  p  +  In  p_) 1 }  *  ♦(In  p)  +  B  (In  p  -  In  p  )2; 

2  A  2  8  1  2 


that  is. 


♦(In  p)  *  1(^  +  ♦  )  -  B  (In  p  -  In  p  )2.  (C6) 

2  1  2  8  1  2 

Subtracting  Eq.  (C5)  from  Eq.  (C4),  on  the  other  hand,  we  find 

A  *  (♦j  -  ♦jJ/dn  px  -  In  p2).  (C7 ) 

Bqs.  (C3) ,  (C6 ) ,  and  (C7)  define  all  the  parameters  for  evaluation  of  ♦(In  p) 
in  Bq.  (Cl),  where  p2  <  p  <  p  .  Once  ♦(In  p)  is  found  for  levels 
bordering  (above  and  below)  the  layer  for  which  T  is  desired,  the  local 
hydrostatic  form 


T  -  d2  -  ♦l)/[R  In  (Pj/p^J 


is  used  to  calculate  layer  temperature. 


B.  Calculation  of  Wind  Components  and  Specific  Humidity 

A  linaar  in  In  p  assumption  is  usad  to  obtain  a  valua  Kin  p)  at  a  dasirad 
prassure  p  from  values  of  X^,  ,  Z^  at  pressure  levels  p^  p2 , 

....  pL,  where  Z  -  u,  v.  In  q.  Assume  than  that 


Z(ln  p)  -  ZCln  p)  +  A(ln  p  -  In  p) 


whara  A  is  a  constant  in  p2  <  p  <  p^  and  In  p  »  “(In  p^  +  In  p2).  Since  Z 
is  linaar  in  In  p, 

az/3  In  p  >  A  >  constant  -  (Z^  -  Z2)/(ln  p1  -  In  p2). 

Than 


and 


Hence 


Z,  -  Z(ln  p)  +  Alin  pt  -  l(ln  p,  +  In  p,)] 

2 


Z2  »  Z(ln  p)  +  Alin  P2  -  l(ln  px  ♦  In  P2>J. 


and 


Z(ln  p)  -  1(Z.  ♦  Z.) 

2  * 


Z(ln  p)  -  1(Z  +  Z  )  + 

2 


*1  -  X2 

_  l In  p  -  l(ln  p  +  In  p,)l. 

In  p.  -  In  p,  2 


This  last  equation  is  the  form  usad  to  interpolate  between  known  values  of  ZJt 
Z2  for  P2  <  P  <  Pj  to  obtain  Z  at  p,  for  Z  *  u,  v,  In  q. 


APPENDIX  D.  Method  of  Obtaining  Satellite  Level  Temperatures 
Prom  Satellite  Layer  Temperatures 


Given  Kg~l  values  of  satellite  layer  temperature  T^  between  layer  interface 
pressures  pfc,  k  *  1,  Kg,  the  Flattery  algorithm  (see  Appendix  A)  is  used  to  get 
preliminary  values  of  the  satellite  level  temperatures  Tfc,  k=*l ,  Kg .  Then  one 
of  these  values  (called  T^at  P^)  is  picked  to  anchor  the  recalculation  of  the 
other  values  in  the  following  way.  Assume  that  the  layer  average  temperature 
is  defined  in  terms  of  the  integral  expression 


"■k  -  ,t  a* 

*x+l 


(Dl) 


where  T  ^  a  +  b  ln  p,  a  linear  function  of  In  p,  and  =  T^.  If  we  denote 

Ap^  *  Pk  -  Pfc+i*  k  »  1,  K  -  1,  then  substituting  the  linear  expression  into 
the  integral  and  performing  the  integration  yields 

\  m  *  *  a>"VIPk(ln  »k  -  11  -  Pktl<ln  "ktl  ‘  in'  lD2) 

A  A 

If  we  choose  to  solve  upward  for  Tfc+1  using  T^  and  Tfe,  then  use 

Tfc  »  a  +  b  In  pk  .  (D3) 

Subtracting  Eq.  (D3)  from  Eq.  (D2)  and  solving  for  b  yields 

b  -  (*k  -  V{Ipk(ln  pk  -  l)  -  pk+i(ln  Pk+1  -  1)1/apk  -  ln  pk)_1  (D4) 

and  from  Eq.  (D3), 

a  »  Tfc  -  b  ln  p^.  (D5) 

Then  T^+1  =  a  +  b  ln  Pk+j  using  these  values  of  a,b. 

If,  on  the  other  hand,  we  wish  to  solve  for  T^  when  T^+k,  T^  are  known 


(solving  downward),  we  use 


Subtracting  Eq.  (D6 )  from  Eq.  (D2)  and  solving  for  b  gives 


b  ■  (Ik  -  W  (It,k<ln  Vl)  -  >Wln  Vi  -  11 -  ln  Vi*'1  <D,) 

and  solving  Eq.  (D6)  for  a  gives 

*  ■  Vi  - b  ln  "k+i  •  (D8> 

Then  T^  *  a  +  b  ln  p^  using  these  values  of  a,  b.  are  the  anchor 

levels  from  which  values  of  T,  at  levels  above  and  below  can  be  evaluated  in  a 

k 

stepwise  fashion. 

This  method  serves  as  an  alternative  to  using  the  Flattery  method  to 
estimate  the  level  temperatures  T^  at  all  levels  from  the  layer  temperatures  T^. 
A  test  was  performed  using  a  radiosonde  sounding  to  simulate  satellite  data  to 
determine  which  of  the  two  methods  produced  the  better  estimates  of  level 
height  from  layer  temperature.  The  following  heights  and  temperatures  were 
taken  from  a  radiosonde  sounding,  for  pressure  levels  usually  included  in  a 
satellite  sounding: 


p  (mb) 

z  W 

T  Ck) 

1000 

100 

291.36 

850 

1479 

284.86 

700 

3080 

275.56 

500 

5716 

259.06 

400 

7364 

246.16 

300 

9389 

232.46 

200 

12023 

212.46 

100 

16235 

206.86 

70 

18376 

207.86 

50 

20423 

209 . 96 

30 

23628 

218.26 

10 

30418 

220.86 

Using  the  above  Z  values,  T  for  each  layer  was  computed  from  the  thickness 
hydrostatically.  These  "satellite  layer  temperatures"  were  converted  to 
"satellite  level  temperatures"  using  two  methods:  (1)  Flattery  method 
entirely,  (2)  Flattery  method  to  get  400  mb  T,  then  integrating  the 
definition  of  layer  mean  temperature  upward  and  downward  as  described 
earlier  in  this  Appendix.  Results  are  shown  in  the  following  table. 


v  (mb) 

Actual  T  C K) 

Method  1  ( °K) 

Method  2  CK) 

1000 

291.36 

293.39 

293.39 

850 

284.86 

286.06 

285.78 

700 

275.56 

276.55 

276.97 

500 

259.06 

258.09 

256.83 

400 

246 . 16 

247.09 

247.09 

300 

232.46 

232.16 

232.86 

200 

212.46 

214.62 

209.11 

100 

206 . 86 

204.57 

205.38 

70 

207.86 

205.88 

204.43 

50 

209.96 

210.85 

211.35 

30 

218.26 

213.94 

217.58 

10 

220.86 

208.67 

201.56 

A  first  guess  surface  pressure  of 

1015.75  mb  at  the 

i  radiosonde  location  was 

used  to  form  the 

basis  for  an  interpolation  of  the 

Method  1  and  2  results  to 

o  layer  pressures 

As  a  reference 

for  comparison. 

the  PT0SIG  routine  is 

used  with  radiosonde  heights  and  temperatures  to  interpolate  Z  quadrat ically 
to  a  levels,  beginning  with  the  first  a  level  above  the  surface  (939.56 
mb)  since  no  extrapolation  is  allowed.  Then  this  o  level  height  is  used  as 
a  base  for  Methods  1  and  2  in  the  computation  of  Z  at  other  levels. 
Temperatures  at  o  layers  from  the  quadratically  interpolated  rawinsonde 
heights  were  also  computed  for  comparison.  All  values  are  shown  in  the 
following  table. 


129 


Sigma  Layer  Temperatures  and  Sigma  Level  Heights 


Ay 


PTOSIG 

METHOD  1 

METHOD  2 

A 

A 

A 

A 

£ 

£ 

Z  T 

Z  T 

Z  T 

939.56 

632.61 

632.61 

632.61 

875.53- 

287.40 

287.39 

287.17 

812.60 

1857.79 

1854.74 

1853.79 

735.47 

278.71 

278.97 

279.22 

660.23 

3549.86 

3551.43 

3551.96 

582.87 

266.41 

266.51 

266.01 

507.87 

5597.21 

5599.50 

5596.20 

443.30 

251.50 

252.16 

251.58 

380.90 

7716.47 

7724.31 

7716.09 

342.31 

239.52 

239.01 

239.38 

304.72 

9282.00 

9286.49 

9280.72 

279.05 

228.20 

229.03 

228.62 

253.94 

10500.67 

10509.61 

10501.62 

228 . 21 

218.06 

220.33 

216.84 

203.15 

11925.91 

11949.70 

11918.89 

177.32 

209.32 

212.88 

208.46 

152.36 

13689.76 

13743.50 

13675.49 

126.36 

206.42 

207.97 

206.64 

101.57 

16141.32 

16213.39 

16129.60 

75.15 

206 . 19 

205.62 

204.62 

50.79 

20327.50 

20388.03 

20284.00 

21.07 

211.99 

212.25 

212.43 

10.00 

30418.00 

30449.84 

30395.41 

In  this 

:  case,  the  Z  values  from  Method  2  are  in  closer  agreement  with  the 

PTOSIG 

(method  of  obtaining  estimates  of  rawinsonde  Z,  T  on  levels,  layers) 

values 

than  are  the  Method  1 

values  at  the  507.87 

mb  level  and  above.  At  the 

lowest 

two  levels  they  are  only  slightly  worse.  On  this  basis,  Method  2 

was 

chosen 

to  process  the  satellite  data  to  produce  a 

level  Z  and  a  layer  T 

estimates,  so  that  their  values  would  be  most  consistent  with  those  for 

rawinsondes . 

APPENDIX  E.  Conversion  of  Cloud  Amount  to  Relative  Humidity 

1 .  AFGWC  Method 

Percent  cloud  amounts  are  converted  to  condensation  pressure  spread  (CPS) 
using  the  curves  given  in  Pig.  El  (data  tables  obtained  from  AFGWC/SDDN, 

Offutt  APB,  ME  68113).  CPS  is  then  converted  to  dew  point  depression  (DPD) 
on  a  pressure  surface  P  (mb)  using 

DPD  =  CPS [4 . 9  +  0.93(P/1000)  +  9 . 0(P/1000) 2 ] _1 . 

Finally,  temperatures  from  an  appropriate  source  (observation  or  analysis)  are 
used  to  convert  DPD  to  relative  humidity. 

2 .  ECMWF  Method 

Total  fractional  cloud  cover  is  estimated  from  observed  relative  humidity 
at  layer  k  using 

RH.  -  RHC. 

CLCk  =  (max[(  ^  ,  0])2 

where 

RHCk  +  1  -  2ok  +  2oZ  +VT  <>k(l  -  3ok  +  2 o£)  . 

°k  =  Pk/P*  • 

Curves  of  relative  humidity  vs.  fractional  cloud  cover  generated  from  this 
expression  are  given  for  several  mandatory  pressures  in  Fig.  E2.  Solving  this 
expression  for  RHk  in  terms  of  fractional  cloud  cover  yields 

RHk  =  RHC  +VCLC  (1  -  RHCk)  for  RHk  >  RHCk  <  1. 

This  computation  requires  a  knowledge  of  ok>  which  in  turn  requires  a 
knowledge  of  the  model  terrain  surface  pressure  at  the  geographic  location  in 
question.  Thus,  in  the  processing  of  the  3DNEPH  data  that  do  not  include 
surface  pressure  information,  one  can  go  only  as  far  as  calculating  v CLC  for 
0.01  <  CLC  <  1,  bilinearly  interpolating  these  to  the  desired  regular 


Pig.  B2  BCMWF  RH  -  Cloud  Amount  Conversion  Scheme 


latitude- longitude  positions,  and  then  later  using  surface  pressures  on  the 
same  positions  to  evaluate  RH^. 

3.  Tibaldi  (from  Chu-Parrish)  Method 

As  used  in  the  surface  weather  observation  conversion  of  cloud  amount  to 
relative  humidity  expression 


RHj  a  Hj  -  Aj  cos[(w/100)  •  X  cloud  amount], 

this  method  simply  uses  percentage  cloud  cover  for  four  cloud  level  categories 
to  estimate  relative  humidity.  The  layers  indexed  by  j  and  the  corresponding 
values  of  and  Aj  are  given  by 

j 

4  high  (P4  -  P5> 

3  middle  (p2  —  P4 

2  low  <p2  -  P3 

1  P.B.L.  (p!  -  p2) 

Fig.  S3  shows  curves  of  relative  humidity  vs.  cloud  amount  (in  eighths)  as 
generated  by  this  equation.  The  layer  pressure  interfaces  are  defined  by 

Pi  -  P* 

p2  *  Pi  -  50  mb 

P3  3  P2  ”  (l/3)(p2  -  P5) 

P4  =  P3  -  (l/3)(p2  -  P5) 

PS  =  300  mb  . 

The  surface  pressure  pA  in  the  analysis  is  the  first  guess  surface  pressure 
on  the  model  terrain.  However,  since  no  information  on  surface  pressure  of 
the  3DMBPH  points  is  included  with  the  cloud  analysis,  it  is  only  possible  to 
compute  cos[(«/100)  •  X  cloud  amount)  at  each  level  in  the  3DWSPH 
processing.  These  values  are  then  bilinearly  interpolated  to  the  desired 
positions  for  use  in  the  analysis.  Then  the  first  guess  surface  pressures  can 
be  used  to  prescribe  the  values  of  p^  -  p$  at  each  observation  or  analysis 
point,  and  these  in  turn  will  dictate  which  values  of  A^ ,  to  use  for 
each  of  the  five  levels  of  the  vertically  compacted  3DNBPH  cloud  amount  values 
to  convert  them  to  relative  humidity. 


M. 

A. 

j 

j 

5T55 

0T10 

0.60 

0.15 

0.75 

0.15 

0.80 

0.20 
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Fig.  E3  HMC  (Chu-Parrish)  RH  -  Cloud  Amount  Conversion  Scheme 


APPENDIX  F.  NMC  Subroutine  GE7PS 


This  subroutine  estimates  the  value  of  In  pA,  where  pA  is  the  surface 
pressure,  from  data  on  known  pressure  levels.  The  geopotential  of  the  level 
where  pA  is  desired,  the  surface  topography,  is  given  and  denoted  by  ♦*. 

Assume  we  are  given  ••••  T^,  T^*  . . . ,  T^  at  the  known 

pressure  levels  p, ,  p_ ,  . . . ,  p. .  Also ,  assume  that  ♦  is  a  quadratic 
function  of  In  p  in  an  individual  layer  p^  <  p  <  pkl  where  here  the 
indices  increase  upward  in  the  atmosphere.  Thus,  we  assume  the  form 

♦(In  p)  =  ♦(In  p)  +  A(ln  p  -  In  p)  +  j  (In  p  -  In  p)2  (FI) 

that  satisfies  the  hydrostatic  equation 

a^/a  In  p  3  -RT  =  A  +  B(ln  p  -  In  p) 

so  that 

a2^/3(ln  p)2  =  -ROT/3  In  p)  -  B  .  (F2) 

The  first  step  is  to  find  the  levels  4^,  that  surround  the  surface 

topography  geopotential  ♦*.  Then  letting 

DH  »  In  pkl  -  In  pK  >  0 

A  =  <.+k  l  -  ^J/DH  <  0 

<  0  lapse 

B  =•  -R(T.  .  -  T.  )/DH  *  0  isothermal 

>  0  inversion 

♦o  =  ♦(IIT7)  -  -(♦k_i  +  ♦jj)  -  |(DH)2. 

If  we  let  XA  =  In  pA  ,  XQ  >  In  p,  then  rewriting  Eq.  (Fl)  gives 


Using  the  quadratic  equation  to  solve  for  X*  -  X  yields 

*  o 


x 


* 


A  +  VD 
B 


where  D  5  A 


2 


-  2B(*  -  >  0.  Choosing  the  +  sign  gives 

O  K 


x 


* 


-  X 

o 


k  +  Vd 

B 


Rationalizing  the  numerator,  using  the  definition  of  D,  and  noting  that 
-A  *  | A |  since  A  <  0  in  all  cases,  yields 


X*  -  Xo  -  2<*0  -  +*)/<  |A|  +  VD) . 


Solving  for  X*  and  using  =  exp(XA)  gives  the  surface  pressure. 
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